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