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