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