0024915: Wrong intersection curves between two cylinders
[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) && ((aU21-aUSurf2l) <= theTol2D) && (aVSurf1f <= aV11) && (aV11 <= aVSurf1l) && (aVSurf2f <= aV21) && (aV21 <= aVSurf2l))
1838         {
1839           if(!aWL1FindStatus)
1840           {
1841             Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
1842
1843             AddBoundaryPoint(theQuad1, theQuad2, aWLine1, anEquationCoeffs, theUVSurf1, theUVSurf2,
1844                         theTol2D, aPeriod, aNulValue, anU1, aU21, aV11, aV11Prev, aV21, aV21Prev, isTheReverse, 1.0, isFound1, isFound2);
1845               
1846             if(isFound1 || isFound2)
1847             {
1848               aWL1FindStatus = 1;
1849             }
1850           }
1851
1852           if((aWL1FindStatus != 2) || (aWLine1->NbPnts() >= 1))
1853           {
1854             if(AddPointIntoWL(theQuad1, theQuad2, isTheReverse, gp_Pnt2d(anU1, aV11), gp_Pnt2d(aU21, aV21), aUSurf1f, aUSurf1l, aPeriod, aWLine1->Curve(), theTol2D))
1855             {
1856               if(!aWL1FindStatus)
1857               {
1858                 aWL1FindStatus = 1;
1859               }
1860             }
1861           }
1862         }
1863         else
1864         {
1865           if(aWL1FindStatus == 1)
1866           {
1867             Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
1868
1869             AddBoundaryPoint(theQuad1, theQuad2, aWLine1, anEquationCoeffs, theUVSurf1, theUVSurf2,
1870                         theTol2D, aPeriod, aNulValue, anU1, aU21, aV11, aV11Prev, aV21, aV21Prev, isTheReverse, 1.0, isFound1, isFound2);
1871
1872             if(isFound1 || isFound2)
1873               aWL1FindStatus = 2; //start a new line
1874           }
1875         }
1876       
1877         if(((aUSurf2f-aU22) <= theTol2D) && ((aU22-aUSurf2l) <= theTol2D) && (aVSurf1f <= aV12) && (aV12 <= aVSurf1l)&& (aVSurf2f <= aV22) && (aV22 <= aVSurf2l))
1878         {
1879           if(!aWL2FindStatus)
1880           {
1881             Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
1882
1883             AddBoundaryPoint(theQuad1, theQuad2, aWLine2, anEquationCoeffs, theUVSurf1, theUVSurf2,
1884                         theTol2D, aPeriod, aNulValue, anU1, aU22, aV12, aV12Prev, aV22, aV22Prev, isTheReverse, -1.0, isFound1, isFound2);
1885               
1886             if(isFound1 || isFound2)
1887             {
1888               aWL2FindStatus = 1;
1889             }
1890           }
1891
1892           if((aWL2FindStatus != 2) || (aWLine2->NbPnts() >= 1))
1893           {
1894             if(AddPointIntoWL(theQuad1, theQuad2, isTheReverse, gp_Pnt2d(anU1, aV12), gp_Pnt2d(aU22, aV22), aUSurf1f, aUSurf1l, aPeriod, aWLine2->Curve(), theTol2D))
1895             {
1896               if(!aWL2FindStatus)
1897               {
1898                 aWL2FindStatus = 1;
1899               }
1900             }
1901           }
1902         }
1903         else
1904         {
1905           if(aWL2FindStatus == 1)
1906           {
1907             Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
1908
1909             AddBoundaryPoint(theQuad1, theQuad2, aWLine2, anEquationCoeffs, theUVSurf1, theUVSurf2,
1910                         theTol2D, aPeriod, aNulValue, anU1, aU22, aV12, aV12Prev, aV22, aV22Prev, isTheReverse, -1.0, isFound1, isFound2);
1911
1912             if(isFound1 || isFound2)
1913               aWL2FindStatus = 2; //start a new line
1914           }
1915         }
1916
1917         aV11Prev = aV11;
1918         aV12Prev = aV12;
1919         aV21Prev = aV21;
1920         aV22Prev = aV22;
1921
1922         if((aWL1FindStatus == 2) || (aWL2FindStatus == 2))
1923         {//current lines are filled. Go to the next lines
1924           anUf = anU1;
1925           break;
1926         }
1927
1928         Standard_Real aFact1 = !IsEqual(sin(aU21 - anEquationCoeffs.mFI2), 0.0) ? 
1929                                   anEquationCoeffs.mK1 * sin(anU1 - anEquationCoeffs.mFIV1) +
1930                                   anEquationCoeffs.mL1 * anEquationCoeffs.mB * sin(aU21 - anEquationCoeffs.mPSIV1) *
1931                                   sin(anU1 - anEquationCoeffs.mFI1)/sin(aU21-anEquationCoeffs.mFI2) : 0.0,
1932                       aFact2 = !IsEqual(sin(aU22-anEquationCoeffs.mFI2), 0.0) ?
1933                                   anEquationCoeffs.mK1 * sin(anU1 - anEquationCoeffs.mFIV1) + 
1934                                   anEquationCoeffs.mL1 * anEquationCoeffs.mB * sin(aU22 - anEquationCoeffs.mPSIV1) *
1935                                   sin(anU1 - anEquationCoeffs.mFI1)/sin(aU22-anEquationCoeffs.mFI2) : 0.0;
1936
1937         Standard_Real aDeltaV1 = (aVSurf1l - aVSurf1f)/IntToReal(aNbPoints);
1938
1939         if((aV11 < aVSurf1f) && (aFact1 < 0.0))
1940         {//Make close to aVSurf1f by increasing anU1 (for the 1st line)
1941           aDeltaV1 = Min(aDeltaV1, Abs(aV11-aVSurf1f));
1942         }
1943
1944         if((aV12 < aVSurf1f) && (aFact2 < 0.0))
1945         {//Make close to aVSurf1f by increasing anU1 (for the 2nd line)
1946           aDeltaV1 = Min(aDeltaV1, Abs(aV12-aVSurf1f));
1947         }
1948
1949         if((aV11 > aVSurf1l) && (aFact1 > 0.0))
1950         {//Make close to aVSurf1l by increasing anU1 (for the 1st line)
1951           aDeltaV1 = Min(aDeltaV1, Abs(aV11-aVSurf1l));
1952         }
1953       
1954         if((aV12 > aVSurf1l) && (aFact2 > 0.0))
1955         {//Make close to aVSurf1l by increasing anU1 (for the 1st line)
1956           aDeltaV1 = Min(aDeltaV1, Abs(aV12-aVSurf1l));
1957         }
1958
1959         Standard_Real aDeltaU1L1 = !IsEqual(aFact1,0.0)? Abs(aDeltaV1/aFact1) : aStepMax,
1960                       aDeltaU1L2 = !IsEqual(aFact2,0.0)? Abs(aDeltaV1/aFact2) : aStepMax;
1961       
1962         const Standard_Real aDeltaU1V1 = Min(aDeltaU1L1, aDeltaU1L2);
1963
1964         ///////////////////////////
1965         aFact1 = !IsEqual(sin(aU21-anEquationCoeffs.mFI2), 0.0) ? 
1966                                   anEquationCoeffs.mK2 * sin(anU1 - anEquationCoeffs.mFIV2) + 
1967                                   anEquationCoeffs.mL2 * anEquationCoeffs.mB * sin(aU21 - anEquationCoeffs.mPSIV2) *
1968                                   sin(anU1 - anEquationCoeffs.mFI1)/sin(aU21 - anEquationCoeffs.mFI2) : 0.0;
1969         aFact2 = !IsEqual(sin(aU22-anEquationCoeffs.mFI2), 0.0) ? 
1970                                   anEquationCoeffs.mK2 * sin(anU1 - anEquationCoeffs.mFIV2) + 
1971                                   anEquationCoeffs.mL2 * anEquationCoeffs.mB * sin(aU22 - anEquationCoeffs.mPSIV2) *
1972                                   sin(anU1 - anEquationCoeffs.mFI1)/sin(aU22 - anEquationCoeffs.mFI2) : 0.0;
1973       
1974         Standard_Real aDeltaV2 = (aVSurf2l - aVSurf2f)/IntToReal(aNbPoints);
1975
1976         if((aV21 < aVSurf2f) && (aFact1 < 0.0))
1977         {//Make close to aVSurf2f by increasing anU1 (for the 1st line)
1978           aDeltaV2 = Min(aDeltaV2, Abs(aV21-aVSurf2f));
1979         }
1980
1981         if((aV22 < aVSurf2f) && (aFact2 < 0.0))
1982         {//Make close to aVSurf1f by increasing anU1 (for the 2nd line)
1983           aDeltaV2 = Min(aDeltaV2, Abs(aV22-aVSurf2f));
1984         }
1985
1986         if((aV21 > aVSurf2l) && (aFact1 > 0.0))
1987         {//Make close to aVSurf1l by increasing anU1 (for the 1st line)
1988           aDeltaV2 = Min(aDeltaV2, Abs(aV21-aVSurf2l));
1989         }
1990
1991         if((aV22 > aVSurf2l) && (aFact2 > 0.0))
1992         {//Make close to aVSurf1l by increasing anU1 (for the 1st line)
1993           aDeltaV2 = Min(aDeltaV2, Abs(aV22-aVSurf1l));
1994         }
1995
1996         aDeltaU1L1 = !IsEqual(aFact1,0.0)? Abs(aDeltaV2/aFact1) : aStepMax;
1997         aDeltaU1L2 = !IsEqual(aFact2,0.0)? Abs(aDeltaV2/aFact2) : aStepMax;
1998
1999         const Standard_Real aDeltaU1V2 = Min(aDeltaU1L1, aDeltaU1L2);
2000
2001         Standard_Real aDeltaU1 = Min(aDeltaU1V1, aDeltaU1V2);
2002
2003         if(aDeltaU1 < aStepMin)
2004           aDeltaU1 = aStepMin;
2005       
2006         if(aDeltaU1 > aStepMax)
2007           aDeltaU1 = aStepMax;
2008
2009         anU1 += aDeltaU1;
2010
2011         const Standard_Real aDiff = anU1 - anUl;
2012         if((0.0 < aDiff) && (aDiff < aDeltaU1-Precision::PConfusion()))
2013           anU1 = anUl;
2014
2015         anUf = anU1;
2016
2017         if(aWLine1->NbPnts() != 1)
2018           isAddedIntoWL1 = Standard_False;
2019
2020         if(aWLine2->NbPnts() != 1)
2021           isAddedIntoWL2 = Standard_False;
2022       }
2023
2024       if((aWLine1->NbPnts() == 1) && (!isAddedIntoWL1))
2025       {
2026         isTheEmpty = Standard_False;
2027         Standard_Real u1, v1, u2, v2;
2028         aWLine1->Point(1).Parameters(u1, v1, u2, v2);
2029         IntPatch_Point aP;
2030         aP.SetParameter(u1);
2031         aP.SetParameters(u1, v1, u2, v2);
2032         aP.SetTolerance(theTol3D);
2033         aP.SetValue(aWLine1->Point(1).Value());
2034
2035         theSPnt.Append(aP);
2036       }
2037       else if(aWLine1->NbPnts() > 1)
2038       {
2039         isTheEmpty = Standard_False;
2040         isAddedIntoWL1 = Standard_True;
2041
2042         SeekAdditionalPoints(theQuad1, theQuad2, aWLine1->Curve(), anEquationCoeffs, aNbPoints, aUSurf2f, aUSurf2l, theTol2D, aPeriod, 1.0, isTheReverse);
2043
2044         aWLine1->ComputeVertexParameters(theTol3D);
2045         theSlin.Append(aWLine1);
2046       }
2047       else
2048       {
2049         isAddedIntoWL1 = Standard_False;
2050       }
2051
2052       if((aWLine2->NbPnts() == 1) && (!isAddedIntoWL2))
2053       {
2054         isTheEmpty = Standard_False;
2055         Standard_Real u1, v1, u2, v2;
2056         aWLine2->Point(1).Parameters(u1, v1, u2, v2);
2057         IntPatch_Point aP;
2058         aP.SetParameter(u1);
2059         aP.SetParameters(u1, v1, u2, v2);
2060         aP.SetTolerance(theTol3D);
2061         aP.SetValue(aWLine2->Point(1).Value());
2062
2063         theSPnt.Append(aP);
2064       }
2065       else if(aWLine2->NbPnts() > 1)
2066       {
2067         isTheEmpty = Standard_False;
2068         isAddedIntoWL2 = Standard_True;
2069
2070         SeekAdditionalPoints(theQuad1, theQuad2, aWLine2->Curve(), anEquationCoeffs, aNbPoints, aUSurf2f, aUSurf2l, theTol2D, aPeriod, -1.0, isTheReverse);
2071
2072         aWLine2->ComputeVertexParameters(theTol3D);
2073         theSlin.Append(aWLine2);
2074       }
2075       else
2076       {
2077         isAddedIntoWL2 = Standard_False;
2078       }
2079     }
2080   }
2081
2082   return Standard_True;
2083 }
2084
2085 //=======================================================================
2086 //function : IntCySp
2087 //purpose  : 
2088 //=======================================================================
2089 Standard_Boolean IntCySp(const IntSurf_Quadric& Quad1,
2090                          const IntSurf_Quadric& Quad2,
2091                          const Standard_Real Tol,
2092                          const Standard_Boolean Reversed,
2093                          Standard_Boolean& Empty,
2094                          Standard_Boolean& Multpoint,
2095                          IntPatch_SequenceOfLine& slin,
2096                          IntPatch_SequenceOfPoint& spnt)
2097
2098 {
2099   Standard_Integer i;
2100
2101   IntSurf_TypeTrans trans1,trans2;
2102   IntAna_ResultType typint;
2103   IntPatch_Point ptsol;
2104   gp_Circ cirsol;
2105
2106   gp_Cylinder Cy;
2107   gp_Sphere Sp;
2108
2109   if (!Reversed) {
2110     Cy = Quad1.Cylinder();
2111     Sp = Quad2.Sphere();
2112   }
2113   else {
2114     Cy = Quad2.Cylinder();
2115     Sp = Quad1.Sphere();
2116   }
2117   IntAna_QuadQuadGeo inter(Cy,Sp,Tol);
2118
2119   if (!inter.IsDone()) {return Standard_False;}
2120
2121   typint = inter.TypeInter();
2122   Standard_Integer NbSol = inter.NbSolutions();
2123   Empty = Standard_False;
2124
2125   switch (typint) {
2126
2127   case IntAna_Empty :
2128     {
2129       Empty = Standard_True;
2130     }
2131     break;
2132
2133   case IntAna_Point :
2134     {
2135       gp_Pnt psol(inter.Point(1));
2136       Standard_Real U1,V1,U2,V2;
2137       Quad1.Parameters(psol,U1,V1);
2138       Quad2.Parameters(psol,U2,V2);
2139       ptsol.SetValue(psol,Tol,Standard_True);
2140       ptsol.SetParameters(U1,V1,U2,V2);
2141       spnt.Append(ptsol);
2142     }
2143     break;
2144
2145   case IntAna_Circle:
2146     {
2147       cirsol = inter.Circle(1);
2148       gp_Vec Tgt;
2149       gp_Pnt ptref;
2150       ElCLib::D1(0.,cirsol,ptref,Tgt);
2151
2152       if (NbSol == 1) {
2153         gp_Vec TestCurvature(ptref,Sp.Location());
2154         gp_Vec Normsp,Normcyl;
2155         if (!Reversed) {
2156           Normcyl = Quad1.Normale(ptref);
2157           Normsp  = Quad2.Normale(ptref);
2158         }
2159         else {
2160           Normcyl = Quad2.Normale(ptref);
2161           Normsp  = Quad1.Normale(ptref);
2162         }
2163         
2164         IntSurf_Situation situcyl;
2165         IntSurf_Situation situsp;
2166         
2167         if (Normcyl.Dot(TestCurvature) > 0.) {
2168           situsp = IntSurf_Outside;
2169           if (Normsp.Dot(Normcyl) > 0.) {
2170             situcyl = IntSurf_Inside;
2171           }
2172           else {
2173             situcyl = IntSurf_Outside;
2174           }
2175         }
2176         else {
2177           situsp = IntSurf_Inside;
2178           if (Normsp.Dot(Normcyl) > 0.) {
2179             situcyl = IntSurf_Outside;
2180           }
2181           else {
2182             situcyl = IntSurf_Inside;
2183           }
2184         }
2185         Handle(IntPatch_GLine) glig;
2186         if (!Reversed) {
2187           glig = new IntPatch_GLine(cirsol, Standard_True, situcyl, situsp);
2188         }
2189         else {
2190           glig = new IntPatch_GLine(cirsol, Standard_True, situsp, situcyl);
2191         }
2192         slin.Append(glig);
2193       }
2194       else {
2195         if (Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)) > 0.0) {
2196           trans1 = IntSurf_Out;
2197           trans2 = IntSurf_In;
2198         }
2199         else {
2200           trans1 = IntSurf_In;
2201           trans2 = IntSurf_Out;
2202         }
2203         Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
2204         slin.Append(glig);
2205
2206         cirsol = inter.Circle(2);
2207         ElCLib::D1(0.,cirsol,ptref,Tgt);
2208         Standard_Real qwe = Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
2209         if(qwe> 0.0000001) {
2210           trans1 = IntSurf_Out;
2211           trans2 = IntSurf_In;
2212         }
2213         else if(qwe<-0.0000001) {
2214           trans1 = IntSurf_In;
2215           trans2 = IntSurf_Out;
2216         }
2217         else { 
2218           trans1=trans2=IntSurf_Undecided; 
2219         }
2220         glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
2221         slin.Append(glig);
2222       }
2223     }
2224     break;
2225     
2226   case IntAna_NoGeometricSolution:
2227     {
2228       gp_Pnt psol;
2229       Standard_Real U1,V1,U2,V2;
2230       IntAna_IntQuadQuad anaint(Cy,Sp,Tol);
2231       if (!anaint.IsDone()) {
2232         return Standard_False;
2233       }
2234
2235       if (anaint.NbPnt()==0 && anaint.NbCurve()==0) {
2236         Empty = Standard_True;
2237       }
2238       else {
2239
2240         NbSol = anaint.NbPnt();
2241         for (i = 1; i <= NbSol; i++) {
2242           psol = anaint.Point(i);
2243           Quad1.Parameters(psol,U1,V1);
2244           Quad2.Parameters(psol,U2,V2);
2245           ptsol.SetValue(psol,Tol,Standard_True);
2246           ptsol.SetParameters(U1,V1,U2,V2);
2247           spnt.Append(ptsol);
2248         }
2249         
2250         gp_Pnt ptvalid,ptf,ptl;
2251         gp_Vec tgvalid;
2252         Standard_Real first,last,para;
2253         IntAna_Curve curvsol;
2254         Standard_Boolean tgfound;
2255         Standard_Integer kount;
2256         
2257         NbSol = anaint.NbCurve();
2258         for (i = 1; i <= NbSol; i++) {
2259           curvsol = anaint.Curve(i);
2260           curvsol.Domain(first,last);
2261           ptf = curvsol.Value(first);
2262           ptl = curvsol.Value(last);
2263
2264           para = last;
2265           kount = 1;
2266           tgfound = Standard_False;
2267           
2268           while (!tgfound) {
2269             para = (1.123*first + para)/2.123;
2270             tgfound = curvsol.D1u(para,ptvalid,tgvalid);
2271             if (!tgfound) {
2272               kount ++;
2273               tgfound = kount > 5;
2274             }
2275           }
2276           Handle(IntPatch_ALine) alig;
2277           if (kount <= 5) {
2278             Standard_Real qwe = tgvalid.DotCross(Quad2.Normale(ptvalid),
2279                                                  Quad1.Normale(ptvalid));
2280             if(qwe> 0.00000001) {
2281               trans1 = IntSurf_Out;
2282               trans2 = IntSurf_In;
2283             }
2284             else if(qwe<-0.00000001) {
2285               trans1 = IntSurf_In;
2286               trans2 = IntSurf_Out;
2287             }
2288             else { 
2289               trans1=trans2=IntSurf_Undecided; 
2290             }
2291             alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
2292           }
2293           else {
2294             alig = new IntPatch_ALine(curvsol,Standard_False);
2295           }
2296           Standard_Boolean TempFalse1a = Standard_False;
2297           Standard_Boolean TempFalse2a = Standard_False;
2298
2299           //-- ptf et ptl : points debut et fin de alig 
2300           
2301           ProcessBounds(alig,slin,Quad1,Quad2,TempFalse1a,ptf,first,
2302                         TempFalse2a,ptl,last,Multpoint,Tol);
2303           slin.Append(alig);
2304         } //-- boucle sur les lignes 
2305       } //-- solution avec au moins une lihne 
2306     }
2307     break;
2308
2309   default:
2310     {
2311       return Standard_False;
2312     }
2313   }
2314   return Standard_True;
2315 }
2316 //=======================================================================
2317 //function : IntCyCo
2318 //purpose  : 
2319 //=======================================================================
2320 Standard_Boolean IntCyCo(const IntSurf_Quadric& Quad1,
2321                          const IntSurf_Quadric& Quad2,
2322                          const Standard_Real Tol,
2323                          const Standard_Boolean Reversed,
2324                          Standard_Boolean& Empty,
2325                          Standard_Boolean& Multpoint,
2326                          IntPatch_SequenceOfLine& slin,
2327                          IntPatch_SequenceOfPoint& spnt)
2328
2329 {
2330   IntPatch_Point ptsol;
2331
2332   Standard_Integer i;
2333
2334   IntSurf_TypeTrans trans1,trans2;
2335   IntAna_ResultType typint;
2336   gp_Circ cirsol;
2337
2338   gp_Cylinder Cy;
2339   gp_Cone     Co;
2340
2341   if (!Reversed) {
2342     Cy = Quad1.Cylinder();
2343     Co = Quad2.Cone();
2344   }
2345   else {
2346     Cy = Quad2.Cylinder();
2347     Co = Quad1.Cone();
2348   }
2349   IntAna_QuadQuadGeo inter(Cy,Co,Tol);
2350
2351   if (!inter.IsDone()) {return Standard_False;}
2352
2353   typint = inter.TypeInter();
2354   Standard_Integer NbSol = inter.NbSolutions();
2355   Empty = Standard_False;
2356
2357   switch (typint) {
2358
2359   case IntAna_Empty : {
2360     Empty = Standard_True;
2361   }
2362     break;
2363
2364   case IntAna_Point :{
2365     gp_Pnt psol(inter.Point(1));
2366     Standard_Real U1,V1,U2,V2;
2367     Quad1.Parameters(psol,U1,V1);
2368     Quad1.Parameters(psol,U2,V2);
2369     ptsol.SetValue(psol,Tol,Standard_True);
2370     ptsol.SetParameters(U1,V1,U2,V2);
2371     spnt.Append(ptsol);
2372   }
2373     break;
2374     
2375   case IntAna_Circle:  {
2376     gp_Vec Tgt;
2377     gp_Pnt ptref;
2378     Standard_Integer j;
2379     Standard_Real qwe;
2380     //
2381     for(j=1; j<=2; ++j) { 
2382       cirsol = inter.Circle(j);
2383       ElCLib::D1(0.,cirsol,ptref,Tgt);
2384       qwe = Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
2385       if(qwe> 0.00000001) {
2386         trans1 = IntSurf_Out;
2387         trans2 = IntSurf_In;
2388       }
2389       else if(qwe<-0.00000001) {
2390         trans1 = IntSurf_In;
2391         trans2 = IntSurf_Out;
2392       }
2393       else { 
2394         trans1=trans2=IntSurf_Undecided; 
2395       }
2396       Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
2397       slin.Append(glig);
2398     }
2399   }
2400     break;
2401     
2402   case IntAna_NoGeometricSolution:    {
2403     gp_Pnt psol;
2404     Standard_Real U1,V1,U2,V2;
2405     IntAna_IntQuadQuad anaint(Cy,Co,Tol);
2406     if (!anaint.IsDone()) {
2407       return Standard_False;
2408     }
2409     
2410     if (anaint.NbPnt() == 0 && anaint.NbCurve() == 0) {
2411       Empty = Standard_True;
2412     }
2413     else {
2414       NbSol = anaint.NbPnt();
2415       for (i = 1; i <= NbSol; i++) {
2416         psol = anaint.Point(i);
2417         Quad1.Parameters(psol,U1,V1);
2418         Quad2.Parameters(psol,U2,V2);
2419         ptsol.SetValue(psol,Tol,Standard_True);
2420         ptsol.SetParameters(U1,V1,U2,V2);
2421         spnt.Append(ptsol);
2422       }
2423       
2424       gp_Pnt ptvalid, ptf, ptl;
2425       gp_Vec tgvalid;
2426       //
2427       Standard_Real first,last,para;
2428       Standard_Boolean tgfound,firstp,lastp,kept;
2429       Standard_Integer kount;
2430       //
2431       //
2432       //IntAna_Curve curvsol;
2433       IntAna_Curve aC;
2434       IntAna_ListOfCurve aLC;
2435       IntAna_ListIteratorOfListOfCurve aIt;
2436       
2437       //
2438       NbSol = anaint.NbCurve();
2439       for (i = 1; i <= NbSol; ++i) {
2440         kept = Standard_False;
2441         //curvsol = anaint.Curve(i);
2442         aC=anaint.Curve(i);
2443         aLC.Clear();
2444         ExploreCurve(Cy, Co, aC, 10.*Tol, aLC);
2445         //
2446         aIt.Initialize(aLC);
2447         for (; aIt.More(); aIt.Next()) {
2448           IntAna_Curve& curvsol=aIt.Value();
2449           //
2450           curvsol.Domain(first, last);
2451           firstp = !curvsol.IsFirstOpen();
2452           lastp  = !curvsol.IsLastOpen();
2453           if (firstp) {
2454             ptf = curvsol.Value(first);
2455           }
2456           if (lastp) {
2457             ptl = curvsol.Value(last);
2458           }
2459           para = last;
2460           kount = 1;
2461           tgfound = Standard_False;
2462           
2463           while (!tgfound) {
2464             para = (1.123*first + para)/2.123;
2465             tgfound = curvsol.D1u(para,ptvalid,tgvalid);
2466             if (!tgfound) {
2467               kount ++;
2468               tgfound = kount > 5;
2469             }
2470           }
2471           Handle(IntPatch_ALine) alig;
2472           if (kount <= 5) {
2473             Standard_Real qwe = tgvalid.DotCross(Quad2.Normale(ptvalid),
2474                                                  Quad1.Normale(ptvalid));
2475             if(qwe> 0.00000001) {
2476               trans1 = IntSurf_Out;
2477               trans2 = IntSurf_In;
2478             }
2479             else if(qwe<-0.00000001) {
2480               trans1 = IntSurf_In;
2481               trans2 = IntSurf_Out;
2482             }
2483             else { 
2484               trans1=trans2=IntSurf_Undecided; 
2485             }
2486             alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
2487             kept = Standard_True;
2488           }
2489           else {
2490             ptvalid = curvsol.Value(para);
2491             alig = new IntPatch_ALine(curvsol,Standard_False);
2492             kept = Standard_True;
2493             //-- cout << "Transition indeterminee" << endl;
2494           }
2495           if (kept) {
2496             Standard_Boolean Nfirstp = !firstp;
2497             Standard_Boolean Nlastp  = !lastp;
2498             ProcessBounds(alig,slin,Quad1,Quad2,Nfirstp,ptf,first,
2499                           Nlastp,ptl,last,Multpoint,Tol);
2500             slin.Append(alig);
2501           }
2502         } // for (; aIt.More(); aIt.Next())
2503       } // for (i = 1; i <= NbSol; ++i) 
2504     }
2505   }
2506     break;
2507     
2508   default:
2509     return Standard_False;
2510     
2511   } // switch (typint)
2512   
2513   return Standard_True;
2514 }
2515 //=======================================================================
2516 //function : ExploreCurve
2517 //purpose  : 
2518 //=======================================================================
2519 Standard_Boolean ExploreCurve(const gp_Cylinder& ,//aCy,
2520                               const gp_Cone& aCo,
2521                               IntAna_Curve& aC,
2522                               const Standard_Real aTol,
2523                               IntAna_ListOfCurve& aLC)
2524                               
2525 {
2526   Standard_Boolean bFind=Standard_False;
2527   Standard_Real aTheta, aT1, aT2, aDst;
2528   gp_Pnt aPapx, aPx;
2529   //
2530   //aC.Dump();
2531   //
2532   aLC.Clear();
2533   aLC.Append(aC);
2534   //
2535   aPapx=aCo.Apex();
2536   //
2537   aC.Domain(aT1, aT2);
2538   //
2539   aPx=aC.Value(aT1);
2540   aDst=aPx.Distance(aPapx);
2541   if (aDst<aTol) {
2542     return bFind;
2543   }
2544   aPx=aC.Value(aT2);
2545   aDst=aPx.Distance(aPapx);
2546   if (aDst<aTol) {
2547     return bFind;
2548   }
2549   //
2550   bFind=aC.FindParameter(aPapx, aTheta);
2551   if (!bFind){
2552     return bFind;
2553   }
2554   //
2555   aPx=aC.Value(aTheta);
2556   aDst=aPx.Distance(aPapx);
2557   if (aDst>aTol) {
2558     return !bFind;
2559   }
2560   //
2561   // need to be splitted at aTheta
2562   IntAna_Curve aC1, aC2;
2563   //
2564   aC1=aC;
2565   aC1.SetDomain(aT1, aTheta);
2566   aC2=aC;
2567   aC2.SetDomain(aTheta, aT2);
2568   //
2569   aLC.Clear();
2570   aLC.Append(aC1);
2571   aLC.Append(aC2);
2572   //
2573   return bFind;
2574 }
2575 //=======================================================================
2576 //function : IsToReverse
2577 //purpose  : 
2578 //=======================================================================
2579 Standard_Boolean IsToReverse(const gp_Cylinder& Cy1,
2580                              const gp_Cylinder& Cy2,
2581                              const Standard_Real Tol)
2582 {
2583   Standard_Boolean bRet;
2584   Standard_Real aR1,aR2, dR, aSc1, aSc2;
2585   //
2586   bRet=Standard_False;
2587   //
2588   aR1=Cy1.Radius();
2589   aR2=Cy2.Radius();
2590   dR=aR1-aR2;
2591   if (dR<0.) {
2592     dR=-dR;
2593   }
2594   if (dR>Tol) {
2595     bRet=aR1>aR2;
2596   }
2597   else {
2598     gp_Dir aDZ(0.,0.,1.);
2599     //
2600     const gp_Dir& aDir1=Cy1.Axis().Direction();
2601     aSc1=aDir1*aDZ;
2602     if (aSc1<0.) {
2603       aSc1=-aSc1;
2604     }
2605     //
2606     const gp_Dir& aDir2=Cy2.Axis().Direction();
2607     aSc2=aDir2*aDZ;
2608     if (aSc2<0.) {
2609       aSc2=-aSc2;
2610     }
2611     //
2612     if (aSc2<aSc1) {
2613       bRet=!bRet;
2614     }
2615   }//if (dR<Tol) {
2616   return bRet;
2617 }