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