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