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