0026937: Eliminate NO_CXX_EXCEPTION macro support
[occt.git] / src / GccAna / GccAna_Circ2d3Tan.cxx
1 // Copyright (c) 1995-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
3 //
4 // This file is part of Open CASCADE Technology software library.
5 //
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
11 //
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
14
15
16 #include <ElCLib.hxx>
17 #include <GccAna_Circ2d3Tan.hxx>
18 #include <GccAna_Circ2dBisec.hxx>
19 #include <GccEnt_BadQualifier.hxx>
20 #include <GccEnt_QualifiedCirc.hxx>
21 #include <GccEnt_QualifiedLin.hxx>
22 #include <GccInt_BCirc.hxx>
23 #include <GccInt_BElips.hxx>
24 #include <GccInt_BHyper.hxx>
25 #include <GccInt_BLine.hxx>
26 #include <GccInt_IType.hxx>
27 #include <gp_Circ2d.hxx>
28 #include <gp_Dir2d.hxx>
29 #include <gp_Lin2d.hxx>
30 #include <gp_Pnt2d.hxx>
31 #include <IntAna2d_AnaIntersection.hxx>
32 #include <IntAna2d_Conic.hxx>
33 #include <IntAna2d_IntPoint.hxx>
34 #include <math_DirectPolynomialRoots.hxx>
35 #include <Standard_OutOfRange.hxx>
36 #include <StdFail_NotDone.hxx>
37 #include <TColStd_Array1OfReal.hxx>
38
39 //=========================================================================
40 //   Creation of a circle tangent to three circles.                        +
41 //=========================================================================
42 GccAna_Circ2d3Tan::
43   GccAna_Circ2d3Tan (const GccEnt_QualifiedCirc& Qualified1,
44                      const GccEnt_QualifiedCirc& Qualified2,
45                      const GccEnt_QualifiedCirc& Qualified3,
46                      const Standard_Real         Tolerance ):
47
48 //=========================================================================
49 //   Initialization of fields.                                           +
50 //=========================================================================
51
52   cirsol(1,16)   ,
53   qualifier1(1,16),
54   qualifier2(1,16),
55   qualifier3(1,16),
56   TheSame1(1,16) ,
57   TheSame2(1,16) ,
58   TheSame3(1,16) ,
59   pnttg1sol(1,16),
60   pnttg2sol(1,16),
61   pnttg3sol(1,16),
62   par1sol(1,16)  ,
63   par2sol(1,16)  ,
64   par3sol(1,16)  ,
65   pararg1(1,16)  ,
66   pararg2(1,16)  ,
67   pararg3(1,16)  
68 {
69
70   gp_Dir2d dirx(1.0,0.0);
71   Standard_Real Tol = Abs(Tolerance);
72   WellDone = Standard_False;
73   NbrSol = 0;
74   if (!(Qualified1.IsEnclosed() || Qualified1.IsEnclosing() || 
75         Qualified1.IsOutside() || Qualified1.IsUnqualified()) ||
76       !(Qualified2.IsEnclosed() || Qualified2.IsEnclosing() || 
77         Qualified2.IsOutside() || Qualified2.IsUnqualified()) ||
78       !(Qualified3.IsEnclosed() || Qualified3.IsEnclosing() || 
79         Qualified3.IsOutside() || Qualified3.IsUnqualified())) {
80     throw GccEnt_BadQualifier();
81       return;
82     }
83   
84 //=========================================================================
85 //   Processing.                                                          +
86 //=========================================================================
87
88   gp_Circ2d Cir1 = Qualified1.Qualified();
89   gp_Circ2d Cir2 = Qualified2.Qualified();
90   gp_Circ2d Cir3 = Qualified3.Qualified();
91   Standard_Real R1 = Cir1.Radius();
92   Standard_Real R2 = Cir2.Radius();
93   Standard_Real R3 = Cir3.Radius();
94   gp_Pnt2d center1(Cir1.Location());
95   gp_Pnt2d center2(Cir2.Location());
96   gp_Pnt2d center3(Cir3.Location());
97   
98   Standard_Real X1 = center1.X();
99   Standard_Real X2 = center2.X();
100   Standard_Real X3 = center3.X();
101   
102   Standard_Real Y1 = center1.Y();
103   Standard_Real Y2 = center2.Y();
104   Standard_Real Y3 = center3.Y();
105
106   gp_XY dir2 = center1.XY() - center2.XY();
107   gp_XY dir3 = center1.XY() - center3.XY();
108
109 //////////
110   if ((Abs(R1 - R2) <= Tolerance && center1.IsEqual(center2, Tolerance)) ||
111       (Abs(R1 - R3) <= Tolerance && center1.IsEqual(center3, Tolerance)) ||
112       (Abs(R2 - R3) <= Tolerance && center2.IsEqual(center3, Tolerance)))
113     return;
114   else {
115     if (Abs(dir2^dir3) <= Tolerance) {
116       Standard_Real Dist1 = center1.Distance(center2);
117       Standard_Real Dist2 = center1.Distance(center3);
118       Standard_Real Dist3 = center2.Distance(center3);
119       if (Abs(Abs(R1 - R2) - Dist1) <= Tolerance) {
120         if (Abs(Abs(R1 - R3) - Dist2) <= Tolerance) {
121           if (Abs(Abs(R2 - R3) - Dist3) <= Tolerance)
122             return;
123         } else if (Abs(R1 + R3 - Dist2) <= Tolerance) {
124           if (Abs(R2 + R3 - Dist3) <= Tolerance)
125             return;
126         }
127       } else if (Abs(R1 + R2 - Dist1) <= Tolerance) {
128         if (Abs(Abs(R1 - R3) - Dist2) <= Tolerance &&
129             Abs(R2 + R3 - Dist3) <= Tolerance) {
130         } else {
131           if (Abs(Abs(R2 - R3) - Dist3) <= Tolerance &&
132               Abs(R1 + R3 - Dist2) <= Tolerance)
133             return;
134         }
135       }
136     }
137   }
138 /////////
139   TColStd_Array1OfReal A2(1, 8), B2(1, 8), C2(1, 8), D2(1, 8), E2(1, 8), F2(1, 8); 
140   TColStd_Array1OfReal A3(1, 8), B3(1, 8), C3(1, 8), D3(1, 8), E3(1, 8), F3(1, 8); 
141   TColStd_Array1OfReal Beta2(1, 8), Gamma2(1, 8), Delta2(1, 8);
142   TColStd_Array1OfReal Beta3(1, 8), Gamma3(1, 8), Delta3(1, 8);
143   Standard_Real a2, b2, c2, d2, e2, f2;
144   Standard_Real a3, b3, c3, d3, e3, f3;
145   Standard_Real A, B, C, D, E;
146   Standard_Boolean IsSame;
147   Standard_Boolean IsTouch;
148   Standard_Integer FirstIndex;
149
150   Standard_Integer i, j, k, l;
151   TColStd_Array1OfReal xSol(1, 64);
152   TColStd_Array1OfReal ySol(1, 64);
153   TColStd_Array1OfReal rSol(1, 16);
154   TColStd_Array1OfInteger FirstSol(1, 9);
155   TColStd_Array1OfReal xSol1(1, 32);
156   TColStd_Array1OfReal ySol1(1, 32);
157   TColStd_Array1OfReal rSol1(1, 32);
158   TColStd_Array1OfInteger FirstSol1(1, 9);
159   Standard_Real x, y, r;
160   Standard_Real m, n, t, s, v;
161   Standard_Real p, q;
162   Standard_Real Epsilon;
163
164   Standard_Integer CurSol;
165
166 //*********************************************************************************************
167 //*********************************************************************************************
168
169 //   Actually we have to find solutions of eight systems of equations:
170 //         _                                          _
171 //        | (X - X1)2 + (Y - Y1)2 = (R - R1)2        | (X - X1)2 + (Y - Y1)2 = (R + R1)2
172 //   1)  <  (X - X2)2 + (Y - Y2)2 = (R - R2)2   2)  <  (X - X2)2 + (Y - Y2)2 = (R - R2)2
173 //        \_(X - X3)2 + (Y - Y3)2 = (R - R3)2        \_(X - X3)2 + (Y - Y3)2 = (R - R3)2
174 //         _                                          _
175 //        | (X - X1)2 + (Y - Y1)2 = (R - R1)2        | (X - X1)2 + (Y - Y1)2 = (R - R1)2
176 //   3)  <  (X - X2)2 + (Y - Y2)2 = (R + R2)2   4)  <  (X - X2)2 + (Y - Y2)2 = (R - R2)2
177 //        \_(X - X3)2 + (Y - Y3)2 = (R - R3)2        \_(X - X3)2 + (Y - Y3)2 = (R + R3)2
178 //         _                                          _
179 //        | (X - X1)2 + (Y - Y1)2 = (R + R1)2        | (X - X1)2 + (Y - Y1)2 = (R + R1)2
180 //   5)  <  (X - X2)2 + (Y - Y2)2 = (R + R2)2   6)  <  (X - X2)2 + (Y - Y2)2 = (R - R2)2
181 //        \_(X - X3)2 + (Y - Y3)2 = (R - R3)2        \_(X - X3)2 + (Y - Y3)2 = (R + R3)2
182 //         _                                          _
183 //        | (X - X1)2 + (Y - Y1)2 = (R - R1)2        | (X - X1)2 + (Y - Y1)2 = (R + R1)2
184 //   7)  <  (X - X2)2 + (Y - Y2)2 = (R + R2)2   8)  <  (X - X2)2 + (Y - Y2)2 = (R + R2)2
185 //        \_(X - X3)2 + (Y - Y3)2 = (R + R3)2        \_(X - X3)2 + (Y - Y3)2 = (R + R3)2
186
187 //   each equation (X - Xi)2 + (Y - Yi)2 = (R +- Ri)2 means that the circle (X,Y,R) is tangent 
188 //   to the circle (Xi,Yi,Ri).
189
190 //   The number of each system is very important. Further index i shows the numer of the system
191
192 //   Further Beta, Gamma and Delta are coefficients of the equation:
193 //                R +- Ri = Beta*X + Gamma*Y + Delta  where i=2 or i=3
194
195 //*********************************************************************************************
196 //*********************************************************************************************
197
198 //   Verification do two circles touch each other or not
199 //   if at least one circle touches other one IsTouch become Standard_Standard_True
200
201   if (Abs((X1 - X2)*(X1 - X2) + (Y1 - Y2)*(Y1 - Y2) - (R1 - R2)*(R1 - R2)) <= Tolerance || 
202       Abs((X1 - X2)*(X1 - X2) + (Y1 - Y2)*(Y1 - Y2) - (R1 + R2)*(R1 + R2)) <= Tolerance || 
203       Abs((X1 - X3)*(X1 - X3) + (Y1 - Y3)*(Y1 - Y3) - (R1 - R3)*(R1 - R3)) <= Tolerance || 
204       Abs((X1 - X3)*(X1 - X3) + (Y1 - Y3)*(Y1 - Y3) - (R1 + R3)*(R1 + R3)) <= Tolerance || 
205       Abs((X2 - X3)*(X2 - X3) + (Y2 - Y3)*(Y2 - Y3) - (R2 - R3)*(R2 - R3)) <= Tolerance || 
206       Abs((X2 - X3)*(X2 - X3) + (Y2 - Y3)*(Y2 - Y3) - (R2 + R3)*(R2 + R3)) <= Tolerance)
207     IsTouch = Standard_True;
208   else
209     IsTouch = Standard_False;
210
211 //   First step:
212 //     We are searching for Beta, Gamma and Delta coefficients
213 //     and also coefficients of the system of second order equations:
214 //     _
215 //    |  a2*x*x +2*b2*x*y + c2*y*y +2*d2*x + 2*e2*y + f2 = 0
216 //   <
217 //    \_ a3*x*x +2*b3*x*y + c3*y*y +2*d3*x + 2*e3*y + f3 = 0   ,
218
219 //     obtained by exclusion of R from source systems.
220
221   for (i = 1; i <= 8; i++) {
222
223 //    _
224 //   | (X - X1)2 + (Y - Y1)2 = (R +- R1)2
225 //  <
226 //   \_(X - X2)2 + (Y - Y2)2 = (R +- R2)2
227
228     if (i == 1 || i == 4 || i == 5 || i == 8) {
229       if (Abs(R1 - R2) > Tolerance) {
230         Beta2(i) = (X1 - X2)/(R1 - R2);
231         Gamma2(i) = (Y1 - Y2)/(R1 - R2);
232         Delta2(i) = (X2*X2 - X1*X1 + Y2*Y2 - Y1*Y1 + (R1 - R2)*(R1 - R2))/(2*(R1 - R2));
233       }
234     } else {
235       Beta2(i) = (X1 - X2)/(R1 + R2);
236       Gamma2(i) = (Y1 - Y2)/(R1 + R2);
237       Delta2(i) = (X2*X2 - X1*X1 + Y2*Y2 - Y1*Y1 + (R1 + R2)*(R1 + R2))/(2*(R1 + R2));
238     }
239     if ((i == 1 || i == 4 || i == 5 || i == 8) &&
240         (Abs(R1 - R2) <= Tolerance)) {
241 //  If R1 = R2
242       A2(i) = 0.;
243       B2(i) = 0.;
244       C2(i) = 0.;
245       D2(i) = X2 - X1;
246       E2(i) = Y2 - Y1;
247       F2(i) = X1*X1 - X2*X2 + Y1*Y1 - Y2*Y2;
248     } else {
249       A2(i) = Beta2(i)*Beta2(i) - 1.;
250       B2(i) = Beta2(i)*Gamma2(i);
251       C2(i) = Gamma2(i)*Gamma2(i) - 1.;
252       D2(i) = Beta2(i)*Delta2(i) + X2;
253       E2(i) = Gamma2(i)*Delta2(i) + Y2;
254       F2(i) = Delta2(i)*Delta2(i) - X2*X2 - Y2*Y2;
255     }
256     
257 //    _
258 //   | (X - X1)2 + (Y - Y1)2 = (R +- R1)2
259 //  <
260 //   \_(X - X3)2 + (Y - Y3)2 = (R +- R3)2
261     
262     if (i == 1 || i == 3 || i == 6 || i == 8) {
263       if (Abs(R1 - R3) > Tolerance) {
264         Beta3(i) = (X1 - X3)/(R1 - R3);
265         Gamma3(i) = (Y1 - Y3)/(R1 - R3);
266         Delta3(i) = (X3*X3 - X1*X1 + Y3*Y3 - Y1*Y1 + (R1 - R3)*(R1 - R3))/(2*(R1 - R3));
267       }
268     } else {
269       Beta3(i) = (X1 - X3)/(R1 + R3);
270       Gamma3(i) = (Y1 - Y3)/(R1 + R3);
271       Delta3(i) = (X3*X3 - X1*X1 + Y3*Y3 - Y1*Y1 + (R1 + R3)*(R1 + R3))/(2*(R1 + R3));
272     }
273     if ((i == 1 || i == 3 || i == 6 || i == 8) &&
274         (Abs(R1 - R3) <= Tolerance)) {
275       A3(i) = 0.;
276       B3(i) = 0.;
277       C3(i) = 0.;
278       D3(i) = X3 - X1;
279       E3(i) = Y3 - Y1;
280       F3(i) = X1*X1 - X3*X3 + Y1*Y1 - Y3*Y3;
281     } else {
282       A3(i) = Beta3(i)*Beta3(i) - 1.;
283       B3(i) = Beta3(i)*Gamma3(i);
284       C3(i) = Gamma3(i)*Gamma3(i) - 1.;
285       D3(i) = Beta3(i)*Delta3(i) + X3;
286       E3(i) = Gamma3(i)*Delta3(i) + Y3;
287       F3(i) = Delta3(i)*Delta3(i) - X3*X3 - Y3*Y3;
288     }
289   }
290
291 //   Second step:
292 //     We are searching for the couple (X,Y) as a solution  of the system:
293 //     _
294 //    |  a2*x*x +2*b2*x*y + c2*y*y +2*d2*x + 2*e2*y + f2 = 0
295 //   <
296 //    \_ a3*x*x +2*b3*x*y + c3*y*y +2*d3*x + 2*e3*y + f3 = 0
297
298   CurSol = 1;
299   for (i = 1; i <= 8; i++) {
300     a2 = A2(i);    a3 = A3(i);
301     b2 = B2(i);    b3 = B3(i);
302     c2 = C2(i);    c3 = C3(i);
303     d2 = D2(i);    d3 = D3(i);
304     e2 = E2(i);    e3 = E3(i);
305     f2 = F2(i);    f3 = F3(i);
306
307     FirstSol(i) = CurSol;
308
309 //  In some cases we know that some systems have no solution in any case due to qualifiers
310     if (((i == 2 || i == 5 || i == 6 || i == 8) && 
311          (Qualified1.IsEnclosed() || Qualified1.IsEnclosing())) ||
312         ((i == 1 || i == 3 || i == 4 || i == 7) && Qualified1.IsOutside()))
313       continue;
314
315     if (((i == 3 || i == 5 || i == 7 || i == 8) && 
316          (Qualified2.IsEnclosed() || Qualified2.IsEnclosing())) ||
317         ((i == 1 || i == 2 || i == 4 || i == 6) && Qualified2.IsOutside()))
318       continue;
319
320     if (((i == 4 || i == 6 || i == 7 || i == 8) && 
321          (Qualified3.IsEnclosed() || Qualified3.IsEnclosing())) ||
322         ((i == 1 || i == 2 || i == 3 || i == 5) && Qualified3.IsOutside()))
323       continue;
324
325 // Check is Cir1 a solution of this system or not
326 // In that case equations are equal to each other
327     if (Abs(a2 - a3) <= Tolerance && Abs(b2 - b3) <= Tolerance && Abs(c2 - c3) <= Tolerance &&
328         Abs(d2 - d3) <= Tolerance && Abs(e2 - e3) <= Tolerance && Abs(f2 - f3) <= Tolerance) {
329       xSol(CurSol) = X1;
330       ySol(CurSol) = Y1;
331       CurSol++;
332       continue;
333     }
334 // 1) a2 = 0
335     if (Abs(a2) <= Tolerance) {
336
337 // 1.1) b2y + d2 = 0
338 //   Searching for solution of the equation Ay2 + By + C = 0
339       A = c2; B = 2.*e2; C = f2;
340       math_DirectPolynomialRoots yRoots(A, B, C);
341       if (yRoots.IsDone() && !yRoots.InfiniteRoots())
342         for (k = 1; k <= yRoots.NbSolutions(); k++) {
343 //   for each y solution:
344           y = yRoots.Value(k);
345 //   Searching for solution of the equation Ax2 + Bx + C = 0
346           if (!(k == 2 && Abs(y - yRoots.Value(1)) <= 10*Tolerance) && 
347               Abs(b2*y + d2) <= b2*Tolerance) {
348             A = a3; B = 2*(b3*y + d3); C = c3*(y*y) + 2*e3*y + f3;
349             math_DirectPolynomialRoots xRoots(A, B, C);
350             if (xRoots.IsDone() && !xRoots.InfiniteRoots())
351               for (j = 1; j <= xRoots.NbSolutions(); j++) {
352                 x = xRoots.Value(j);
353                 if (!(j == 2 && Abs(x - xRoots.Value(1)) <= 10*Tolerance)) {
354                   xSol(CurSol) = x;
355                   ySol(CurSol) = y;
356                   CurSol++;
357                 }
358               }
359           }
360         }
361
362 // 1.2) b2y + d2 != 0
363       A = a3*c2*c2 - 4*b2*(b3*c2 - b2*c3);
364       B = 4*a3*c2*e2 - 4*b3*(c2*d2 + 2*b2*e2) + 4*b2*(2*c3*d2 - c2*d3 + 2*b2*e3);
365       C = 2*a3*(c2*f2 + 2*e2*e2) - 4*b3*(b2*f2 + 2*e2*d2) + 4*c3*d2*d2 - 4*d3*(c2*d2 + 2*b2*e2) 
366         + 16*b2*e3*d2 + 4*b2*b2*f3;
367       D = 4*a3*e2*f2 - 4*b3*d2*f2 - 4*d3*(b2*f2 + 2*d2*e2) + 8*d2*d2*e3 + 8*b2*d2*f3;
368       E = a3*f2*f2 - 4*d2*d3*f2 + 4*d2*d2*f3;
369
370 //   Searching for solution of the equation Ay4 + By3 + Cy2 + Dy + E = 0
371 // Special case: one circle touches other 
372       if (IsTouch) {
373 // Derivation of the equation Ay4 + By3 + Cy2 + Dy + E
374         math_DirectPolynomialRoots yRoots1(4*A, 3*B, 2*C, D);
375         if (yRoots1.IsDone() && !yRoots1.InfiniteRoots())
376           for (k = 1; k <= yRoots1.NbSolutions(); k++) {
377             y = yRoots1.Value(k);
378 // Check if this value is already catched
379             IsSame = Standard_False;
380             for (l = 1; l < k; l++)
381               if (Abs(y - yRoots1.Value(l)) <= 10*Tolerance) IsSame = Standard_True;
382             
383             Epsilon = (Abs((Abs((Abs(4*A*y) + Abs(3*B))*y) + Abs(2*C))*y) + Abs(D));
384             if (Abs((((A*y + B)*y + C)*y + D)*y + E) <= Epsilon*Tolerance) {
385               if (!IsSame && Abs(b2*y + d2) > b2*Tolerance) {
386                 x = -(c2*(y*y) + 2*e2*y + f2)/(2*(b2*y + d2));
387                 xSol(CurSol) = x;
388                 ySol(CurSol) = y;
389                 CurSol++;
390               }
391             }
392           }
393       }
394
395       math_DirectPolynomialRoots yRoots1(A, B, C, D, E);
396       if (yRoots1.IsDone() && !yRoots1.InfiniteRoots())
397         for (k = 1; k <= yRoots1.NbSolutions(); k++) {
398           y = yRoots1.Value(k);
399 // Check if this value is already catched
400           IsSame = Standard_False;
401           FirstIndex = (i == 1) ? 1 : FirstSol(i);
402           for (l = FirstIndex; l < CurSol; l++)
403             if (Abs(y - ySol(l)) <= 10*Tolerance) IsSame = Standard_True;
404
405           if (!IsSame && Abs(b2*y + d2) > b2*Tolerance) {
406             x = -(c2*(y*y) + 2*e2*y + f2)/(2*(b2*y + d2));
407             xSol(CurSol) = x;
408             ySol(CurSol) = y;
409             CurSol++;
410           }
411         }
412     } else {
413 // 2) a2 != 0
414 // Coefficients of the equation     (sy + v)Sqrt(p2 - q) = (my2 + ny + t)
415       m = 2*a3*b2*b2/(a2*a2) - 2*b2*b3/a2 - a3*c2/a2 + c3;
416       n = 4*a3*b2*d2/(a2*a2) - 2*b3*d2/a2 - 2*b2*d3/a2 - 2*a3*e2/a2 + 2*e3;
417       t = 2*a3*d2*d2/(a2*a2) - 2*d2*d3/a2 - a3*f2/a2 + f3;
418       s = 2*b3 - 2*a3*b2/a2;
419       v = 2*d3 - 2*d2*a3/a2;
420
421 //------------------------------------------
422 // If s = v = 0
423       if (Abs(s) <= Tolerance && Abs(v) <= Tolerance) {
424         math_DirectPolynomialRoots yRoots(m, n, t);
425         if (yRoots.IsDone() && !yRoots.InfiniteRoots())
426           for (k = 1; k <= yRoots.NbSolutions(); k++) {
427 //   for each y solution:
428             y = yRoots.Value(k);
429
430             p = -(b2*y + d2)/a2;
431             q = (c2*(y*y) + 2*e2*y + f2)/a2;
432             Epsilon = 2.*(Abs((b2*b2 + Abs(a2*c2))*y) + Abs(b2*d2) + Abs(a2*e2))/(a2*a2);
433             if (!(k == 2 && Abs(y - yRoots.Value(1)) <= 10*Tolerance) &&
434                 p*p - q >= -Epsilon*Tolerance) {
435               A = a2; 
436               B = 2*(b2*y + d2);
437               C = c2*y*y + 2*e2*y + f2;
438               math_DirectPolynomialRoots xRoots(A, B, C);
439               if (xRoots.IsDone() && !xRoots.InfiniteRoots())
440                 for (l = 1; l <= xRoots.NbSolutions(); l++) {
441 //   for each x solution:
442                   x = xRoots.Value(l);
443
444                   if (!(l == 2 && Abs(x - xRoots.Value(1)) <= 10*Tolerance)) {
445                     xSol(CurSol) = x;
446                     ySol(CurSol) = y;
447                     CurSol++;
448                   }
449                 }
450             }
451           }
452       } else {
453 //------------------------------------------
454 // If (s*y + v) != 0
455
456         A = s*s*(b2*b2 - a2*c2) - m*m*a2*a2;
457         B = 2*s*v*(b2*b2 - a2*c2) + 2*s*s*(b2*d2 - a2*e2) - 2*m*n*a2*a2;
458         C = v*v*(b2*b2 - a2*c2) + 4*s*v*(b2*d2 - a2*e2) + s*s*(d2*d2 - a2*f2) - a2*a2*(2*m*t + n*n);
459         D = 2*v*v*(b2*d2 - a2*e2) + 2*s*v*(d2*d2 - a2*f2) - 2*n*t*a2*a2;
460         E = v*v*(d2*d2 - a2*f2) - t*t*a2*a2;
461
462 //   Searching for solution of the equation Ay4 + By3 + Cy2 + Dy + E = 0
463 // Special case: one circle touches other 
464         if (IsTouch) {
465 // Derivation of the equation Ay4 + By3 + Cy2 + Dy + E
466           math_DirectPolynomialRoots yRoots1(4*A, 3*B, 2*C, D);
467           if (yRoots1.IsDone() && !yRoots1.InfiniteRoots())
468             for (k = 1; k <= yRoots1.NbSolutions(); k++) {
469               y = yRoots1.Value(k);
470
471               p = -(b2*y + d2)/a2;
472               q = (c2*(y*y) + 2*e2*y + f2)/a2;
473
474 // Check if this value is already catched
475               IsSame = Standard_False;
476               FirstIndex = (i == 1) ? 1 : FirstSol(i);
477               for (l = FirstIndex; l < CurSol; l++)
478                 if (Abs(y - ySol(l)) <= 10*Tolerance) IsSame = Standard_True;
479
480               Epsilon = (Abs((Abs((Abs(4*A*y) + Abs(3*B))*y) + Abs(2*C))*y) + Abs(D));
481               if (Abs((((A*y + B)*y + C)*y + D)*y + E) <= Epsilon*Tolerance) {
482
483                 Epsilon = 2.*(Abs((b2*b2 + Abs(a2*c2))*y) + Abs(b2*d2) + Abs(a2*e2))/(a2*a2);
484                 if (!IsSame && p*p - q >= -Epsilon*Tolerance) {
485                   A = a2; 
486                   B = 2*(b2*y + d2);
487                   C = c2*y*y + 2*e2*y + f2;
488                   math_DirectPolynomialRoots xRoots(A, B, C);
489                   if (xRoots.IsDone() && !xRoots.InfiniteRoots())
490                     for (l = 1; l <= xRoots.NbSolutions(); l++) {
491 //   for each x solution:
492                       x = xRoots.Value(l);
493
494                       if (!(l == 2 && Abs(x - xRoots.Value(1)) <= 10*Tolerance)) {
495                         xSol(CurSol) = x;
496                         ySol(CurSol) = y;
497                         CurSol++;
498                       }
499                     }
500                 }
501               }
502             }
503         }
504
505         math_DirectPolynomialRoots yRoots(A, B, C, D, E);
506         if (yRoots.IsDone() && !yRoots.InfiniteRoots())
507           for (k = 1; k <= yRoots.NbSolutions(); k++) {
508 //   for each y solution:
509             y = yRoots.Value(k);
510
511             p = -(b2*y + d2)/a2;
512             q = (c2*(y*y) + 2*e2*y + f2)/a2;
513
514 // Check if this value is already catched
515             IsSame = Standard_False;
516             for (l = 1; l < k; l++)
517               if (Abs(y - yRoots.Value(l)) <= 10*Tolerance) IsSame = Standard_True;
518
519             Epsilon = 2.*(Abs((b2*b2 + Abs(a2*c2))*y) + Abs(b2*d2) + Abs(a2*e2))/(a2*a2);
520             if (!IsSame && p*p - q >= -Epsilon*Tolerance) {
521               A = a2; 
522               B = 2*(b2*y + d2);
523               C = c2*y*y + 2*e2*y + f2;
524               math_DirectPolynomialRoots xRoots(A, B, C);
525               if (xRoots.IsDone() && !xRoots.InfiniteRoots())
526                 for (l = 1; l <= xRoots.NbSolutions(); l++) {
527 //   for each x solution:
528                   x = xRoots.Value(l);
529
530                   if (!(l == 2 && Abs(x - xRoots.Value(1)) <= 10*Tolerance)) {
531                     xSol(CurSol) = x;
532                     ySol(CurSol) = y;
533                     CurSol++;
534                   }
535                 }
536             }
537           }
538       }
539     }
540   }
541   FirstSol(9) = CurSol;
542
543 //  Third step:
544 //    Check of couples (X,Y) and searching for R. R must be great than 0
545   CurSol = 1;
546   for (i = 1; i <= 8; i++) {
547     FirstSol1(i) = CurSol;
548     for (j = FirstSol(i); j < FirstSol(i + 1); j++) {
549       x = xSol(j);
550       y = ySol(j);
551 // in some cases when R1 = R2 :
552       if ((i == 1 || i == 4 || i == 5 || i == 8) && (Abs(R1 - R2) <= Tolerance)) {
553         if (i == 1 || i == 4) {
554           r = R1 + Sqrt((x - X1)*(x - X1) + (y - Y1)*(y - Y1));
555           Epsilon = 10*(2*Abs(r - R2) + Abs(x - X2) + Abs(y - Y2));
556           if (Abs((r - R2)*(r - R2) - (x - X2)*(x - X2) - (y - Y2)*(y - Y2)) <= 
557               Epsilon*Tolerance) {
558             xSol1(CurSol) = x;
559             ySol1(CurSol) = y;
560             rSol1(CurSol) = r;
561             CurSol++;
562           }
563           r = R1 - Sqrt((x - X1)*(x - X1) + (y - Y1)*(y - Y1));
564           Epsilon = 10*(2*Abs(r - R2) + Abs(x - X2) + Abs(y - Y2));
565           if ((r > Tolerance) && 
566               (Abs((r - R2)*(r - R2) - (x - X2)*(x - X2) - (y - Y2)*(y - Y2)) <= 
567                Epsilon*Tolerance)) {
568             xSol1(CurSol) = x;
569             ySol1(CurSol) = y;
570             rSol1(CurSol) = r;
571             CurSol++;
572           }
573         } else {
574 //      i == 5 || i == 8
575           r =  - R1 + Sqrt((x - X1)*(x - X1) + (y - Y1)*(y - Y1));
576           if (r > Tolerance) {
577             xSol1(CurSol) = x;
578             ySol1(CurSol) = y;
579             rSol1(CurSol) = r;
580             CurSol++;
581           }
582         }
583       } else {
584 // Other cases
585         if (i == 1 || i == 4) {
586           r = R2 + Beta2(i)*x + Gamma2(i)*y + Delta2(i);
587           if (r > Tolerance) {
588             xSol1(CurSol) = x;
589             ySol1(CurSol) = y;
590             rSol1(CurSol) = r;
591             CurSol++;
592           }
593         }
594         if (i == 5 || i == 8) {
595           r = -R2 - Beta2(i)*x - Gamma2(i)*y - Delta2(i);
596           if (r > Tolerance) {
597             xSol1(CurSol) = x;
598             ySol1(CurSol) = y;
599             rSol1(CurSol) = r;
600             CurSol++;
601           }
602         }
603         if (i == 3 || i == 7) {
604           r = - R2 + Beta2(i)*x + Gamma2(i)*y + Delta2(i);
605           if (r > Tolerance) {
606             xSol1(CurSol) = x;
607             ySol1(CurSol) = y;
608             rSol1(CurSol) = r;
609             CurSol++;
610           }
611         }
612         if (i == 2 || i == 6) {
613           r = R2 - Beta2(i)*x - Gamma2(i)*y - Delta2(i);
614           if (r > Tolerance) {
615             xSol1(CurSol) = x;
616             ySol1(CurSol) = y;
617             rSol1(CurSol) = r;
618             CurSol++;
619           }
620         }
621       }
622     }
623   }
624   FirstSol1(9) = CurSol;
625 //   Fourth step
626 //     Check of triplets (X,Y,R).
627   CurSol = 1;
628   for (i = 1; i <= 8; i++) {
629     FirstSol(i) = CurSol;
630     for (j = FirstSol1(i); j < FirstSol1(i + 1); j++) {
631       x = xSol1(j);
632       y = ySol1(j);
633       r = rSol1(j);
634 // in some cases when R1 = R3 :
635       if ((i == 1 || i == 3 || i == 6 || i == 8) && Abs(R1 - R3) <= Tolerance) {
636         if (i == 1 || i == 3) {
637           Epsilon = 10*(2*Abs(r - R3) + Abs(x - X3) + Abs(y - Y3));
638           if (Abs((r - R3)*(r - R3) - (x - X3)*(x - X3) - (y - Y3)*(y - Y3)) <= 
639               Epsilon*Tolerance) {
640             xSol(CurSol) = x;
641             ySol(CurSol) = y;
642             rSol(CurSol) = r;
643             CurSol++;
644           }
645         } else {
646 //      i == 6 || i == 8
647           Epsilon = 10*(2*(r + R3) + Abs(x - X3) + Abs(y - Y3));
648           if (Abs((r + R3)*(r + R3) - (x - X3)*(x - X3) - (y - Y3)*(y - Y3)) <= 
649               Epsilon*Tolerance) {
650             xSol(CurSol) = x;
651             ySol(CurSol) = y;
652             rSol(CurSol) = r;
653             CurSol++;
654           }
655         }
656       } else {
657 // Other cases
658         Epsilon = 10*(Abs(Beta3(i)) + Abs(Gamma3(i)) + 1.);
659         if (i == 1 || i == 3)
660           if (Abs(R3 + Beta3(i)*x + Gamma3(i)*y + Delta3(i) - r) <= Epsilon*Tolerance) {
661             xSol(CurSol) = x;
662             ySol(CurSol) = y;
663             rSol(CurSol) = r;
664             CurSol++;
665           }
666         if (i == 6 || i == 8)
667           if (Abs(R3 + Beta3(i)*x + Gamma3(i)*y + Delta3(i) + r) <= Epsilon*Tolerance) {
668             xSol(CurSol) = x;
669             ySol(CurSol) = y;
670             rSol(CurSol) = r;
671             CurSol++;
672           }
673         if (i == 4 || i == 7)
674           if (Abs(Beta3(i)*x + Gamma3(i)*y + Delta3(i) - r - R3) <= Epsilon*Tolerance) {
675             xSol(CurSol) = x;
676             ySol(CurSol) = y;
677             rSol(CurSol) = r;
678             CurSol++;
679           }
680         if (i == 2 || i == 5)
681           if (Abs(r - R3 + Beta3(i)*x + Gamma3(i)*y + Delta3(i)) <= Epsilon*Tolerance) {
682             xSol(CurSol) = x;
683             ySol(CurSol) = y;
684             rSol(CurSol) = r;
685             CurSol++;
686           }
687       }
688     }
689   }
690   FirstSol(9) = CurSol;
691
692 //  Fifth step:
693 //    We have found all solutions. We have to calculate some parameters for each one.
694   for (i = 1 ; i <= 8; i++) {
695     for (j = FirstSol(i); j < FirstSol(i + 1); j++) {
696
697       if ((Qualified1.IsEnclosed() && rSol(j) > R1) ||
698           (Qualified1.IsEnclosing() && rSol(j) < R1))
699         continue;
700       if ((Qualified2.IsEnclosed() && rSol(j) > R2) ||
701           (Qualified2.IsEnclosing() && rSol(j) < R2))
702         continue;
703       if ((Qualified3.IsEnclosed() && rSol(j) > R3) ||
704           (Qualified3.IsEnclosing() && rSol(j) < R3))
705         continue;
706
707       NbrSol++;
708
709 // RLE, avoid out of range
710       if (NbrSol > cirsol.Upper()) NbrSol = cirsol.Upper();
711       
712       gp_Pnt2d Center = gp_Pnt2d(xSol(j), ySol(j));
713
714       cirsol(NbrSol) = gp_Circ2d(gp_Ax2d(Center,dirx),rSol(j));
715
716 //   ==========================================================
717       Standard_Real distcc1 = Center.Distance(center1);
718       if (!Qualified1.IsUnqualified())
719         qualifier1(NbrSol) = Qualified1.Qualifier();
720       else if (Abs(distcc1 + rSol(j) - R1) <= Tol)
721         qualifier1(NbrSol) = GccEnt_enclosed;
722       else if (Abs(distcc1 - R1 - rSol(j)) <= Tol)
723         qualifier1(NbrSol) = GccEnt_outside;
724       else qualifier1(NbrSol) = GccEnt_enclosing;
725
726       Standard_Real distcc2 = Center.Distance(center1);
727       if (!Qualified2.IsUnqualified())
728         qualifier2(NbrSol) = Qualified2.Qualifier();
729       else if (Abs(distcc2 + rSol(j) - R2) <= Tol)
730         qualifier2(NbrSol) = GccEnt_enclosed;
731       else if (Abs(distcc2 - R2 - rSol(j)) <= Tol)
732         qualifier2(NbrSol) = GccEnt_outside;
733       else qualifier2(NbrSol) = GccEnt_enclosing;
734
735       Standard_Real distcc3 = Center.Distance(center1);
736       if (!Qualified3.IsUnqualified())
737         qualifier3(NbrSol) = Qualified3.Qualifier();
738       else if (Abs(distcc3 + rSol(j) - R3) <= Tol)
739         qualifier3(NbrSol) = GccEnt_enclosed;
740       else if (Abs(distcc3 - R3 - rSol(j)) <= Tol)
741         qualifier3(NbrSol) = GccEnt_outside;
742       else qualifier3(NbrSol) = GccEnt_enclosing;
743
744 //   ==========================================================
745
746       if (Center.Distance(Cir1.Location()) <= Tolerance)
747         TheSame1(NbrSol) = 1;
748       else {
749         TheSame1(NbrSol) = 0;
750         gp_Dir2d dc;
751         if ((i == 2 || i == 5 || i == 6 || i == 8) || rSol(j) > R1)
752           dc.SetXY(Cir1.Location().XY() - Center.XY());
753         else
754           dc.SetXY(Center.XY() - Cir1.Location().XY());
755         pnttg1sol(NbrSol)=gp_Pnt2d(Center.XY() + rSol(j)*dc.XY());
756         par1sol(NbrSol)=ElCLib::Parameter(cirsol(NbrSol),
757                                           pnttg1sol(NbrSol));
758         pararg1(NbrSol)=ElCLib::Parameter(Cir1,pnttg1sol(NbrSol));
759       }
760
761       if (Center.Distance(Cir2.Location()) <= Tolerance)
762         TheSame2(NbrSol) = 1;
763       else {
764         TheSame2(NbrSol) = 0;
765         gp_Dir2d dc;
766         if ((i == 3 || i == 5 || i == 7 || i == 8) || rSol(j) > R2)
767           dc.SetXY(Cir2.Location().XY() - Center.XY());
768         else
769           dc.SetXY(Center.XY() - Cir2.Location().XY());
770         pnttg2sol(NbrSol)=gp_Pnt2d(Center.XY() + rSol(j)*dc.XY());
771         par2sol(NbrSol)=ElCLib::Parameter(cirsol(NbrSol),
772                                           pnttg2sol(NbrSol));
773         pararg2(NbrSol)=ElCLib::Parameter(Cir2,pnttg2sol(NbrSol));
774       }
775
776       if (Center.Distance(Cir3.Location()) <= Tolerance)
777         TheSame3(NbrSol) = 1;
778       else {
779         TheSame3(NbrSol) = 0;
780         gp_Dir2d dc;
781         if ((i == 4 || i == 6 || i == 7 || i == 8) || rSol(j) > R3)
782           dc.SetXY(Cir3.Location().XY() - Center.XY());
783         else
784           dc.SetXY(Center.XY() - Cir3.Location().XY());
785         pnttg3sol(NbrSol)=gp_Pnt2d(Center.XY() + rSol(j)*dc.XY());
786         par3sol(NbrSol)=ElCLib::Parameter(cirsol(NbrSol),
787                                           pnttg3sol(NbrSol));
788         pararg3(NbrSol)=ElCLib::Parameter(Cir3,pnttg3sol(NbrSol));
789       }
790     }
791   }
792   WellDone = Standard_True;
793 }
794
795 //=========================================================================
796
797 Standard_Boolean GccAna_Circ2d3Tan::
798    IsDone () const {
799    return WellDone;
800  }
801
802 Standard_Integer GccAna_Circ2d3Tan::
803    NbSolutions () const {
804    return NbrSol;
805  }
806
807 gp_Circ2d GccAna_Circ2d3Tan::
808    ThisSolution (const Standard_Integer Index) const 
809 {
810   if (!WellDone)
811     throw StdFail_NotDone();
812   
813   if (Index <= 0 ||Index > NbrSol)
814     throw Standard_OutOfRange();
815   
816   return cirsol(Index); 
817 }
818
819 void GccAna_Circ2d3Tan::
820   WhichQualifier(const Standard_Integer Index   ,
821                        GccEnt_Position& Qualif1 ,
822                        GccEnt_Position& Qualif2 ,
823                        GccEnt_Position& Qualif3 ) const
824 {
825   if (!WellDone) { throw StdFail_NotDone(); }
826    else if (Index <= 0 ||Index > NbrSol) { throw Standard_OutOfRange(); }
827    else {
828      Qualif1 = qualifier1(Index);
829      Qualif2 = qualifier2(Index);
830      Qualif3 = qualifier3(Index);
831    }
832 }
833
834 void GccAna_Circ2d3Tan::
835    Tangency1 (const Standard_Integer Index,
836               Standard_Real& ParSol,
837               Standard_Real& ParArg,
838               gp_Pnt2d& PntSol) const {
839    if (!WellDone) {
840      throw StdFail_NotDone();
841    }
842    else if (Index <= 0 ||Index > NbrSol) {
843      throw Standard_OutOfRange();
844    }
845    else {
846      if (TheSame1(Index) == 0) {
847        ParSol = par1sol(Index);
848        ParArg = pararg1(Index);
849        PntSol = gp_Pnt2d(pnttg1sol(Index));
850      }
851      else { throw StdFail_NotDone(); }
852    }
853  }
854
855 void GccAna_Circ2d3Tan::
856    Tangency2 (const Standard_Integer Index,
857               Standard_Real& ParSol,
858               Standard_Real& ParArg,
859               gp_Pnt2d& PntSol) const{
860    if (!WellDone) {
861      throw StdFail_NotDone();
862    }
863    else if (Index <= 0 ||Index > NbrSol) {
864      throw Standard_OutOfRange();
865    }
866    else {
867      if (TheSame2(Index) == 0) {
868        ParSol = par2sol(Index);
869        ParArg = pararg2(Index);
870        PntSol = gp_Pnt2d(pnttg2sol(Index));
871      }
872      else { throw StdFail_NotDone(); }
873    }
874  }
875
876 void GccAna_Circ2d3Tan::
877    Tangency3 (const Standard_Integer Index,
878               Standard_Real& ParSol,
879               Standard_Real& ParArg,
880               gp_Pnt2d& PntSol) const{
881    if (!WellDone) {
882      throw StdFail_NotDone();
883    }
884    else if (Index <= 0 ||Index > NbrSol) {
885      throw Standard_OutOfRange();
886    }
887    else {
888      if (TheSame3(Index) == 0) {
889        ParSol = par3sol(Index);
890        ParArg = pararg3(Index);
891        PntSol = gp_Pnt2d(pnttg3sol(Index));
892      }
893      else { throw StdFail_NotDone(); }
894    }
895  }
896
897 Standard_Boolean GccAna_Circ2d3Tan::
898    IsTheSame1 (const Standard_Integer Index) const
899 {
900   if (!WellDone)
901     throw StdFail_NotDone();
902   
903   if (Index <= 0 ||Index > NbrSol)
904     throw Standard_OutOfRange();
905   
906   if (TheSame1(Index) == 0)
907     return Standard_False;
908   
909   return Standard_True;
910 }
911
912 Standard_Boolean GccAna_Circ2d3Tan::
913    IsTheSame2 (const Standard_Integer Index) const
914 {
915   if (!WellDone)
916     throw StdFail_NotDone();
917
918   if (Index <= 0 ||Index > NbrSol)
919     throw Standard_OutOfRange();
920
921   if (TheSame2(Index) == 0)
922     return Standard_False;
923
924   return Standard_True;
925 }
926
927 Standard_Boolean GccAna_Circ2d3Tan::
928    IsTheSame3 (const Standard_Integer Index) const
929 {
930   if (!WellDone)
931     throw StdFail_NotDone();
932   
933   if (Index <= 0 ||Index > NbrSol)
934     throw Standard_OutOfRange();
935
936   if (TheSame3(Index) == 0) 
937     return Standard_False;
938   
939   return Standard_True;
940 }
941
942