src/Extrema/Extrema_ExtElC.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 <Extrema_ExtElC.hxx>
  18 #include <Extrema_ExtElC2d.hxx>
  19 #include <Extrema_ExtPElC.hxx>
  20 #include <Extrema_POnCurv.hxx>
  21 #include <gp_Ax1.hxx>
  22 #include <gp_Ax2.hxx>
  23 #include <gp_Ax3.hxx>
  24 #include <gp_Circ.hxx>
  25 #include <gp_Circ2d.hxx>
  26 #include <gp_Dir.hxx>
  27 #include <gp_Elips.hxx>
  28 #include <gp_Hypr.hxx>
  29 #include <gp_Lin.hxx>
  30 #include <gp_Lin2d.hxx>
  31 #include <gp_Parab.hxx>
  32 #include <gp_Pln.hxx>
  33 #include <gp_Pnt.hxx>
  34 #include <IntAna_QuadQuadGeo.hxx>
  35 #include <IntAna2d_AnaIntersection.hxx>
  36 #include <math_DirectPolynomialRoots.hxx>
  37 #include <math_TrigonometricFunctionRoots.hxx>
  38 #include <Precision.hxx>
  39 #include <Standard_OutOfRange.hxx>
  40 #include <StdFail_NotDone.hxx>
  41
  42 #include <stdio.h>
  43 static
  44   void RefineDir(gp_Dir& aDir);
  45
  46 //=======================================================================
  47 //class    : ExtremaExtElC_TrigonometricRoots
  48 //purpose  :
  49 //==  Classe Interne (Donne des racines classees d un polynome trigo)
  50 //==  Code duplique avec IntAna_IntQuadQuad.cxx (lbr le 26 mars 98)
  51 //==  Solution fiable aux problemes de coefficients proches de 0
  52 //==  avec essai de rattrapage si coeff<1.e-10 (jct le 27 avril 98)
  53 //=======================================================================
  54 class ExtremaExtElC_TrigonometricRoots {
  55  private:
  56   Standard_Real Roots[4];
  57   Standard_Boolean done;
  58   Standard_Integer NbRoots;
  59   Standard_Boolean infinite_roots;
  60  public:
  61   ExtremaExtElC_TrigonometricRoots(const Standard_Real CC,
  62                                    const Standard_Real SC,
  63                                    const Standard_Real C,
  64                                    const Standard_Real S,
  65                                    const Standard_Real Cte,
  66                                    const Standard_Real Binf,
  67                                    const Standard_Real Bsup);
  68   //
  69   Standard_Boolean IsDone() {
  70     return done;
  71   }
  72   //
  73   Standard_Boolean IsARoot(Standard_Real u) {
  74     Standard_Real PIpPI, aEps;
  75     //
  76     aEps=RealEpsilon();
  77     PIpPI = M_PI + M_PI;
  78     for(Standard_Integer i=0 ; i<NbRoots; i++) {
  79       if(Abs(u - Roots[i])<=aEps) {
  80         return Standard_True ;
  81       }
  82       if(Abs(u - Roots[i]-PIpPI)<=aEps) {
  83         return Standard_True;
  84       }
  85     }
  86     return Standard_False;
  87   }
  88   //
  89   Standard_Integer NbSolutions() {
  90     if(!done) {
  91       throw StdFail_NotDone();
  92     }
  93     return NbRoots;
  94   }
  95   //
  96   Standard_Boolean InfiniteRoots() {
  97     if(!done) {
  98       throw StdFail_NotDone();
  99     }
 100     return infinite_roots;
 101   }
 102   //
 103   Standard_Real Value(const Standard_Integer& n) {
 104     if((!done)||(n>NbRoots)) {
 105       throw StdFail_NotDone();
 106     }
 107     return Roots[n-1];
 108   }
 109 };
 110 //=======================================================================
 111 //function : ExtremaExtElC_TrigonometricRoots
 112 //purpose  :
 113 //=======================================================================
 114 ExtremaExtElC_TrigonometricRoots::
 115   ExtremaExtElC_TrigonometricRoots(const Standard_Real CC,
 116                                    const Standard_Real SC,
 117                                    const Standard_Real C,
 118                                    const Standard_Real S,
 119                                    const Standard_Real Cte,
 120                                    const Standard_Real Binf,
 121                                    const Standard_Real Bsup)
 122 {
 123   Standard_Integer i, nbessai;
 124   Standard_Real cc ,sc, c, s, cte;
 125   //
 126   nbessai = 1;
 127   cc = CC;
 128   sc = SC;
 129   c = C;
 130   s = S;
 131   cte = Cte;
 132   done=Standard_False;
 133   while (nbessai<=2 && !done) {
 134     //-- F= AA*CN*CN+2*BB*CN*SN+CC*CN+DD*SN+EE;
 135     math_TrigonometricFunctionRoots MTFR(cc,sc,c,s,cte,Binf,Bsup);
 136     //
 137     if(MTFR.IsDone()) {
 138       done=Standard_True;
 139       if(MTFR.InfiniteRoots()) {
 140         infinite_roots=Standard_True;
 141       }
 142       else { //else #1
 143         Standard_Boolean Triee;
 144         Standard_Integer j, SvNbRoots;
 145         Standard_Real aTwoPI, aMaxCoef, aPrecision;
 146         //
 147         aTwoPI=M_PI+M_PI;
 148         NbRoots=MTFR.NbSolutions();
 149         for(i=0;i<NbRoots;++i) {
 150           Roots[i]=MTFR.Value(i+1);
 151           if(Roots[i]<0.) {
 152             Roots[i]=Roots[i]+aTwoPI;
 153           }
 154           if(Roots[i]>aTwoPI) {
 155             Roots[i]=Roots[i]-aTwoPI;
 156           }
 157         }
 158         //
 159         //-- La recherche directe donne n importe quoi.
 160         aMaxCoef = Max(CC,SC);
 161         aMaxCoef = Max(aMaxCoef,C);
 162         aMaxCoef = Max(aMaxCoef,S);
 163         aMaxCoef = Max(aMaxCoef,Cte);
 164         aPrecision = Max(1.e-8, 1.e-12*aMaxCoef);
 165
 166         SvNbRoots=NbRoots;
 167         for(i=0; i<SvNbRoots; ++i) {
 168           Standard_Real y;
 169           Standard_Real co=cos(Roots[i]);
 170           Standard_Real si=sin(Roots[i]);
 171           y=co*(CC*co + (SC+SC)*si + C) + S*si + Cte;
 172           // modified by OCC  Tue Oct  3 18:43:00 2006
 173           if(Abs(y)>aPrecision) {
 174             NbRoots--;
 175             Roots[i]=1000.0;
 176           }
 177         }
 178         //
 179         do {
 180           Standard_Real t;
 181           //
 182           Triee=Standard_True;
 183           for(i=1, j=0; i<SvNbRoots; ++i, ++j) {
 184             if(Roots[i]<Roots[j]) {
 185               Triee=Standard_False;
 186               t=Roots[i];
 187               Roots[i]=Roots[j];
 188               Roots[j]=t;
 189             }
 190           }
 191         }
 192         while(!Triee);
 193         //
 194         infinite_roots=Standard_False;
 195         if(NbRoots==0) { //--!!!!! Detect case Pol = Cte ( 1e-50 ) !!!!
 196           if((Abs(CC) + Abs(SC) + Abs(C) + Abs(S)) < 1e-10) {
 197             if(Abs(Cte) < 1e-10)  {
 198               infinite_roots=Standard_True;
 199             }
 200           }
 201         }
 202       } // else #1
 203     } // if(MTFR.IsDone()) {
 204     else {
 205         // try to set very small coefficients to ZERO
 206       if (Abs(CC)<1e-10) {
 207         cc = 0.0;
 208       }
 209       if (Abs(SC)<1e-10) {
 210         sc = 0.0;
 211       }
 212       if (Abs(C)<1e-10) {
 213         c = 0.0;
 214       }
 215       if (Abs(S)<1e-10){
 216         s = 0.0;
 217       }
 218       if (Abs(Cte)<1e-10){
 219         cte = 0.0;
 220       }
 221       nbessai++;
 222     }
 223   } // while (nbessai<=2 && !done) {
 224 }
 225
 226 //=======================================================================
 227 //function : Extrema_ExtElC
 228 //purpose  :
 229 //=======================================================================
 230 Extrema_ExtElC::Extrema_ExtElC ()
 231 {
 232   myDone = Standard_False;
 233   myIsPar = Standard_False;
 234   myNbExt = 0;
 235
 236   for (Standard_Integer i = 0; i < 6; i++)
 237   {
 238     mySqDist[i] = RealLast();
 239   }
 240 }
 241 //=======================================================================
 242 //function : Extrema_ExtElC
 243 //purpose  :
 244 //=======================================================================
 245 Extrema_ExtElC::Extrema_ExtElC (const gp_Lin& theC1,
 246                                 const gp_Lin& theC2,
 247                                 const Standard_Real)
 248 // Function:
 249 //   Find min distance between 2 straight lines.
 250
 251 // Method:
 252 //   Let D1 and D2, be 2 directions of straight lines C1 and C2.
 253 //   2 cases are considered:
 254 //   1- if Angle(D1,D2) < AngTol, straight lines are parallel.
 255 //      The distance is the distance between a point of C1 and the straight line C2.
 256 //   2- if Angle(D1,D2) > AngTol:
 257 //      Let P1=C1(u1) and P2=C2(u2).
 258 //      Then we must find u1 and u2 such as the distance P1P2 is minimal.
 259 //      Target function is:
 260 //        F(u1, u2) = ((L1x+D1x*u1)-(L2x+D2x*u2))^2 +
 261 //                    ((L1y+D1y*u1)-(L2y+D2y*u2))^2 +
 262 //                    ((L1z+D1z*u1)-(L2z+D2z*u2))^2 --> min,
 263 //      where L1 and L2 are lines locations, D1 and D2 are lines directions.
 264 //    Let simplify the function F
 265
 266 //      F(u1, u2) = (D1x*u1-D2x*u2-Lx)^2 + (D1y*u1-D2y*u2-Ly)^2 + (D1z*u1-D2z*u2-Lz)^2,
 267 //    where L is a vector L1L2.
 268
 269 //    In point of minimum, the condition
 270 //      {dF/du1 = 0
 271 //      {dF/du2 = 0
 272
 273 //    must be satisfied.
 274
 275 //      dF/du1 = 2*D1x*(D1x*u1-D2x*u2-Lx) +
 276 //               2*D1y*(D1y*u1-D2y*u2-Ly) +
 277 //               2*D1z*(D1z*u1-D2z*u2-Lz) =
 278 //             = 2*((D1^2)*u1-(D1.D2)*u2-L.D1) =
 279 //             = 2*(u1-(D1.D2)*u2-L.D1)
 280 //      dF/du2 = -2*D2x*(D1x*u1-D2x*u2-Lx) -
 281 //                2*D2y*(D1y*u1-D2y*u2-Ly) -
 282 //                2*D2z*(D1z*u1-D2z*u2-Lz)=
 283 //             = -2*((D2.D1)*u1-(D2^2)*u2-(D2.L)) =
 284 //             = -2*((D2.D1)*u1-u2-(D2.L))
 285
 286 //   Consequently, we have two-equation system with two variables:
 287
 288 //     {u1 - (D1.D2)*u2 = L.D1
 289 //     {(D1.D2)*u1 - u2 = L.D2
 290
 291 //   This system has one solution if (D1.D2)^2 != 1
 292 //   (if straight lines are not parallel).
 293 {
 294   myDone = Standard_False;
 295   myNbExt = 0;
 296   myIsPar = Standard_False;
 297
 298   const gp_Dir &aD1 = theC1.Position().Direction(),
 299                &aD2 = theC2.Position().Direction();
 300   const Standard_Real aCosA = aD1.Dot(aD2);
 301   const Standard_Real aSqSinA = 1.0 - aCosA * aCosA;
 302   Standard_Real aU1 = 0.0, aU2 = 0.0;
 303   if (aSqSinA < gp::Resolution() || aD1.IsParallel(aD2, Precision::Angular()))
 304   {
 305     myIsPar = Standard_True;
 306   }
 307   else
 308   {
 309     const gp_XYZ aL1L2 = theC2.Location().XYZ() - theC1.Location().XYZ();
 310     const Standard_Real aD1L = aD1.XYZ().Dot(aL1L2),
 311                         aD2L = aD2.XYZ().Dot(aL1L2);
 312     aU1 = (aD1L - aCosA * aD2L) / aSqSinA;
 313     aU2 = (aCosA * aD1L - aD2L) / aSqSinA;
 314
 315     myIsPar = Precision::IsInfinite(aU1) || Precision::IsInfinite(aU2);
 316   }
 317
 318   if (myIsPar)
 319   {
 320     mySqDist[0] = theC2.SquareDistance(theC1.Location());
 321     myNbExt = 1;
 322     myDone = Standard_True;
 323     return;
 324   }
 325
 326   // Here myIsPar == Standard_False;
 327
 328   const gp_Pnt aP1(ElCLib::Value(aU1, theC1)),
 329                aP2(ElCLib::Value(aU2, theC2));
 330   mySqDist[myNbExt] = aP1.SquareDistance(aP2);
 331   myPoint[myNbExt][0] = Extrema_POnCurv(aU1, aP1);
 332   myPoint[myNbExt][1] = Extrema_POnCurv(aU2, aP2);
 333   myNbExt = 1;
 334   myDone = Standard_True;
 335 }
 336
 337 //=======================================================================
 338 //function : PlanarLineCircleExtrema
 339 //purpose  :
 340 //=======================================================================
 341 Standard_Boolean Extrema_ExtElC::PlanarLineCircleExtrema(const gp_Lin& theLin,
 342                                                          const gp_Circ& theCirc)
 343 {
 344   const gp_Dir &aDirC = theCirc.Axis().Direction(),
 345                &aDirL = theLin.Direction();
 346
 347   if (Abs(aDirC.Dot(aDirL)) > Precision::Angular())
 348     return Standard_False;
 349
 350   //The line is in the circle-plane completely
 351   //(or parallel to the circle-plane).
 352   //Therefore, we are looking for extremas and
 353   //intersections in 2D-space.
 354
 355   const gp_XYZ &aCLoc = theCirc.Location().XYZ();
 356   const gp_XYZ &aDCx = theCirc.Position().XDirection().XYZ(),
 357                &aDCy = theCirc.Position().YDirection().XYZ();
 358
 359   const gp_XYZ &aLLoc = theLin.Location().XYZ();
 360   const gp_XYZ &aLDir = theLin.Direction().XYZ();
 361
 362   const gp_XYZ aVecCL(aLLoc - aCLoc);
 363
 364   //Center of 2D-circle
 365   const gp_Pnt2d aPC(0.0, 0.0);
 366
 367   gp_Ax22d aCircAxis(aPC, gp_Dir2d(1.0, 0.0), gp_Dir2d(0.0, 1.0));
 368   gp_Circ2d aCirc2d(aCircAxis, theCirc.Radius());
 369
 370   gp_Pnt2d aPL(aVecCL.Dot(aDCx), aVecCL.Dot(aDCy));
 371   gp_Dir2d aDL(aLDir.Dot(aDCx), aLDir.Dot(aDCy));
 372   gp_Lin2d aLin2d(aPL, aDL);
 373
 374   // Extremas
 375   Extrema_ExtElC2d anExt2d(aLin2d, aCirc2d, Precision::Confusion());
 376   //Intersections
 377   IntAna2d_AnaIntersection anInters(aLin2d, aCirc2d);
 378
 379   myDone = anExt2d.IsDone() || anInters.IsDone();
 380
 381   if (!myDone)
 382     return Standard_True;
 383
 384   const Standard_Integer aNbExtr = anExt2d.NbExt();
 385   const Standard_Integer aNbSol = anInters.NbPoints();
 386
 387   const Standard_Integer aNbSum = aNbExtr + aNbSol;
 388   for (Standard_Integer anExtrID = 1; anExtrID <= aNbSum; anExtrID++)
 389   {
 390     const Standard_Integer aDelta = anExtrID - aNbExtr;
 391
 392     Standard_Real aLinPar = 0.0, aCircPar = 0.0;
 393
 394     if (aDelta < 1)
 395     {
 396       Extrema_POnCurv2d aPLin2d, aPCirc2d;
 397       anExt2d.Points(anExtrID, aPLin2d, aPCirc2d);
 398       aLinPar = aPLin2d.Parameter();
 399       aCircPar = aPCirc2d.Parameter();
 400     }
 401     else
 402     {
 403       aLinPar = anInters.Point(aDelta).ParamOnFirst();
 404       aCircPar = anInters.Point(aDelta).ParamOnSecond();
 405     }
 406
 407     const gp_Pnt aPOnL(ElCLib::LineValue(aLinPar, theLin.Position())),
 408                  aPOnC(ElCLib::CircleValue(aCircPar,
 409                                             theCirc.Position(), theCirc.Radius()));
 410
 411     mySqDist[myNbExt] = aPOnL.SquareDistance(aPOnC);
 412     myPoint[myNbExt][0].SetValues(aLinPar, aPOnL);
 413     myPoint[myNbExt][1].SetValues(aCircPar, aPOnC);
 414     myNbExt++;
 415   }
 416
 417   return Standard_True;
 418 }
 419
 420 //=======================================================================
 421 //function : Extrema_ExtElC
 422 //purpose  :
 423 // Find extreme distances between straight line C1 and circle C2.
 424 //
 425 //Method:
 426 //   Let P1=C1(u1) and P2=C2(u2) be two solution points
 427 //        D the direction of straight line C1
 428 //      T tangent at point P2;
 429 //  Then, ( P1P2.D = 0. (1)
 430 //         ( P1P2.T = 0. (2)
 431 //  Let O1 and O2 be the origins of C1 and C2;
 432 //  Then, (1) <=> (O1P2-u1*D).D = 0.         as O1P1 = u1*D
 433 //           <=> u1 = O1P2.D                as D.D = 1.
 434 //         (2) <=> P1O2.T = 0.                as O2P2.T = 0.
 435 //             <=> ((P2O1.D)D+O1O2).T = 0.    as P1O1 = -u1*D = (P2O1.D)D
 436 //           <=> (((P2O2+O2O1).D)D+O1O2).T = 0.
 437 //           <=> ((P2O2.D)(D.T)+((O2O1.D)D-O2O1).T = 0.
 438 //  We are in the reference of the circle; let:
 439 //         Cos = Cos(u2) and Sin = Sin(u2),
 440 //         P2 (R*Cos,R*Sin,0.),
 441 //         T (-R*Sin,R*Cos,0.),
 442 //       D (Dx,Dy,Dz),
 443 //       V (Vx,Vy,Vz) = (O2O1.D)D-O2O1;
 444 //  Then, the equation by Cos and Sin is as follows:
 445 //    -(2*R*R*Dx*Dy)   * Cos**2  +       A1
 446 //   R*R*(Dx**2-Dy**2) * Cos*Sin +    2* A2
 447 //         R*Vy        * Cos     +       A3
 448 //      -R*Vx        * Sin     +       A4
 449 //      R*R*Dx*Dy                = 0.    A5
 450 //Use the algorithm math_TrigonometricFunctionRoots to solve this equation.
 451 //=======================================================================
 452 Extrema_ExtElC::Extrema_ExtElC (const gp_Lin& C1,
 453                                 const gp_Circ& C2,
 454                                 const Standard_Real)
 455 {
 456   Standard_Real Dx,Dy,Dz,aRO2O1, aTolRO2O1;
 457   Standard_Real R, A1, A2, A3, A4, A5, aTol;
 458   gp_Dir x2, y2, z2, D, D1;
 459   //
 460   myIsPar = Standard_False;
 461   myDone = Standard_False;
 462   myNbExt = 0;
 463
 464   if (PlanarLineCircleExtrema(C1, C2))
 465   {
 466     return;
 467   }
 468
 469   // Calculate T1 in the reference of the circle ...
 470   D = C1.Direction();
 471   D1 = D;
 472   x2 = C2.XAxis().Direction();
 473   y2 = C2.YAxis().Direction();
 474   z2 = C2.Axis().Direction();
 475   Dx = D.Dot(x2);
 476   Dy = D.Dot(y2);
 477   Dz = D.Dot(z2);
 478   //
 479   D.SetCoord(Dx, Dy, Dz);
 480   RefineDir(D);
 481   D.Coord(Dx, Dy, Dz);
 482   //
 483   // Calcul de V dans le repere du cercle:
 484   gp_Pnt O1 = C1.Location();
 485   gp_Pnt O2 = C2.Location();
 486   gp_Vec O2O1 (O2,O1);
 487   //
 488   aTolRO2O1=gp::Resolution();
 489   aRO2O1=O2O1.Magnitude();
 490   if (aRO2O1 > aTolRO2O1) {
 491     gp_Dir aDO2O1;
 492     //
 493     O2O1.Multiply(1./aRO2O1);
 494     aDO2O1.SetCoord(O2O1.Dot(x2), O2O1.Dot(y2), O2O1.Dot(z2));
 495     RefineDir(aDO2O1);
 496     O2O1.SetXYZ(aRO2O1*aDO2O1.XYZ());
 497   }
 498   else {
 499     O2O1.SetCoord(O2O1.Dot(x2), O2O1.Dot(y2), O2O1.Dot(z2));
 500   }
 501   //
 502   gp_XYZ Vxyz = (D.XYZ()*(O2O1.Dot(D)))-O2O1.XYZ();
 503   //
 504   //modified by NIZNHY-PKV Tue Mar 20 10:36:38 2012
 505   /*
 506   R = C2.Radius();
 507   A5 = R*R*Dx*Dy;
 508   A1 = -2.*A5;
 509   A2 = R*R*(Dx*Dx-Dy*Dy)/2.;
 510   A3 = R*Vxyz.Y();
 511   A4 = -R*Vxyz.X();
 512   //
 513   aTol=1.e-12;
 514   //
 515
 516   if(fabs(A5) <= aTol) {
 517     A5 = 0.;
 518   }
 519   if(fabs(A1) <= aTol) {
 520     A1 = 0.;
 521   }
 522   if(fabs(A2) <= aTol) {
 523     A2 = 0.;
 524   }
 525   if(fabs(A3) <= aTol) {
 526     A3 = 0.;
 527   }
 528   if(fabs(A4) <= aTol) {
 529     A4 = 0.;
 530   }
 531   */
 532   //
 533   aTol=1.e-12;
 534   // Calculate the coefficients of the equation by Cos and Sin ...
 535   // [divided by R]
 536   R = C2.Radius();
 537   A5 = R*Dx*Dy;
 538   A1 = -2.*A5;
 539   A2 = 0.5*R*(Dx*Dx-Dy*Dy);// /2.;
 540   A3 = Vxyz.Y();
 541   A4 = -Vxyz.X();
 542   //
 543   if (A1>=-aTol && A1<=aTol) {
 544     A1 = 0.;
 545   }
 546   if (A2>=-aTol && A2<=aTol) {
 547     A2 = 0.;
 548   }
 549   if (A3>=-aTol && A3<=aTol) {
 550     A3 = 0.;
 551   }
 552   if (A4>=-aTol && A4<=aTol) {
 553     A4 = 0.;
 554   }
 555   if (A5>=-aTol && A5<=aTol) {
 556     A5 = 0.;
 557   }
 558   //modified by NIZNHY-PKV Tue Mar 20 10:36:40 2012t
 559   //
 560   ExtremaExtElC_TrigonometricRoots Sol(A1, A2, A3, A4, A5, 0., M_PI+M_PI);
 561   if (!Sol.IsDone()) {
 562     return;
 563   }
 564   if (Sol.InfiniteRoots()) {
 565     myIsPar = Standard_True;
 566     mySqDist[0] = R*R;
 567     myNbExt = 1;
 568     myDone = Standard_True;
 569     return;
 570   }
 571   // Storage of solutions ...
 572   Standard_Integer NoSol, NbSol;
 573   Standard_Real U1,U2;
 574   gp_Pnt P1,P2;
 575   //
 576   NbSol = Sol.NbSolutions();
 577   for (NoSol=1; NoSol<=NbSol; ++NoSol) {
 578     U2 = Sol.Value(NoSol);
 579     P2 = ElCLib::Value(U2,C2);
 580     U1 = (gp_Vec(O1,P2)).Dot(D1);
 581     P1 = ElCLib::Value(U1,C1);
 582     mySqDist[myNbExt] = P1.SquareDistance(P2);
 583     //modified by NIZNHY-PKV Wed Mar 21 08:11:33 2012f
 584     //myPoint[myNbExt][0] = Extrema_POnCurv(U1,P1);
 585     //myPoint[myNbExt][1] = Extrema_POnCurv(U2,P2);
 586     myPoint[myNbExt][0].SetValues(U1,P1);
 587     myPoint[myNbExt][1].SetValues(U2,P2);
 588     //modified by NIZNHY-PKV Wed Mar 21 08:11:36 2012t
 589     myNbExt++;
 590   }
 591   myDone = Standard_True;
 592 }
 593 //=======================================================================
 594 //function : Extrema_ExtElC
 595 //purpose  :
 596 //=======================================================================
 597 Extrema_ExtElC::Extrema_ExtElC (const gp_Lin& C1,
 598                                 const gp_Elips& C2)
 599 {
 600 /*-----------------------------------------------------------------------------
 601 Function:
 602   Find extreme distances between straight line C1 and ellipse C2.
 603
 604 Method:
 605   Let P1=C1(u1) and P2=C2(u2) two solution points
 606         D the direction of straight line C1
 607         T the tangent to point P2;
 608   Then, ( P1P2.D = 0. (1)
 609          ( P1P2.T = 0. (2)
 610   Let O1 and O2 be the origins of C1 and C2;
 611   Then, (1) <=> (O1P2-u1*D).D = 0.        as O1P1 = u1*D
 612              <=> u1 = O1P2.D              as D.D = 1.
 613          (2) <=> P1O2.T = 0.              as O2P2.T = 0.
 614              <=> ((P2O1.D)D+O1O2).T = 0.  as P1O1 = -u1*D = (P2O1.D)D
 615              <=> (((P2O2+O2O1).D)D+O1O2).T = 0.
 616              <=> ((P2O2.D)(D.T)+((O2O1.D)D-O2O1).T = 0.
 617   We are in the reference of the ellipse; let:
 618          Cos = Cos(u2) and Sin = Sin(u2),
 619          P2 (MajR*Cos,MinR*Sin,0.),
 620          T (-MajR*Sin,MinR*Cos,0.),
 621          D (Dx,Dy,Dz),
 622          V (Vx,Vy,Vz) = (O2O1.D)D-O2O1;
 623   Then,  get the following equation by Cos and Sin:
 624     -(2*MajR*MinR*Dx*Dy)             * Cos**2  +
 625    (MajR*MajR*Dx**2-MinR*MinR*Dy**2) * Cos*Sin +
 626          MinR*Vy                     * Cos     +
 627        - MajR*Vx                     * Sin     +
 628       MinR*MajR*Dx*Dy                = 0.
 629   Use algorithm math_TrigonometricFunctionRoots to solve this equation.
 630 -----------------------------------------------------------------------------*/
 631   myIsPar = Standard_False;
 632   myDone = Standard_False;
 633   myNbExt = 0;
 634
 635 // Calculate T1 the reference of the ellipse ...
 636   gp_Dir D = C1.Direction();
 637   gp_Dir D1 = D;
 638   gp_Dir x2, y2, z2;
 639   x2 = C2.XAxis().Direction();
 640   y2 = C2.YAxis().Direction();
 641   z2 = C2.Axis().Direction();
 642   Standard_Real Dx = D.Dot(x2);
 643   Standard_Real Dy = D.Dot(y2);
 644   Standard_Real Dz = D.Dot(z2);
 645   D.SetCoord(Dx,Dy,Dz);
 646
 647 // Calculate V ...
 648   gp_Pnt O1 = C1.Location();
 649   gp_Pnt O2 = C2.Location();
 650   gp_Vec O2O1 (O2,O1);
 651   O2O1.SetCoord(O2O1.Dot(x2), O2O1.Dot(y2), O2O1.Dot(z2));
 652   gp_XYZ Vxyz = (D.XYZ()*(O2O1.Dot(D)))-O2O1.XYZ();
 653
 654 // Calculate the coefficients of the equation by Cos and Sin ...
 655   Standard_Real MajR = C2.MajorRadius();
 656   Standard_Real MinR = C2.MinorRadius();
 657   Standard_Real A5 = MajR*MinR*Dx*Dy;
 658   Standard_Real A1 = -2.*A5;
 659   Standard_Real R2 = MajR*MajR;
 660   Standard_Real r2 = MinR*MinR;
 661   Standard_Real A2 =(R2*Dx*Dx -r2*Dy*Dy -R2 +r2)/2.0;
 662   Standard_Real A3 = MinR*Vxyz.Y();
 663   Standard_Real A4 = -MajR*Vxyz.X();
 664   //
 665   Standard_Real aEps=1.e-12;
 666   //
 667   if(fabs(A5) <= aEps) A5 = 0.;
 668   if(fabs(A1) <= aEps) A1 = 0.;
 669   if(fabs(A2) <= aEps) A2 = 0.;
 670   if(fabs(A3) <= aEps) A3 = 0.;
 671   if(fabs(A4) <= aEps) A4 = 0.;
 672   //
 673   ExtremaExtElC_TrigonometricRoots Sol(A1,A2,A3,A4,A5,0.,M_PI+M_PI);
 674   if (!Sol.IsDone()) { return; }
 675   //
 676   if (Sol.InfiniteRoots()) {
 677     myIsPar = Standard_True;
 678     gp_Pnt aP = ElCLib::EllipseValue(0., C2.Position(), C2.MajorRadius(), C2.MinorRadius());
 679     mySqDist[0] = C1.SquareDistance(aP);
 680     myNbExt = 1;
 681     myDone = Standard_True;
 682     return;
 683   }
 684
 685 // Storage of solutions ...
 686   gp_Pnt P1,P2;
 687   Standard_Real U1,U2;
 688   Standard_Integer NbSol = Sol.NbSolutions();
 689   for (Standard_Integer NoSol = 1; NoSol <= NbSol; NoSol++) {
 690     U2 = Sol.Value(NoSol);
 691     P2 = ElCLib::Value(U2,C2);
 692     U1 = (gp_Vec(O1,P2)).Dot(D1);
 693     P1 = ElCLib::Value(U1,C1);
 694     mySqDist[myNbExt] = P1.SquareDistance(P2);
 695     myPoint[myNbExt][0] = Extrema_POnCurv(U1,P1);
 696     myPoint[myNbExt][1] = Extrema_POnCurv(U2,P2);
 697     myNbExt++;
 698   }
 699   myDone = Standard_True;
 700 }
 701
 702 //=======================================================================
 703 //function : Extrema_ExtElC
 704 //purpose  :
 705 //=======================================================================
 706 Extrema_ExtElC::Extrema_ExtElC (const gp_Lin& C1,
 707                                 const gp_Hypr& C2)
 708 {
 709 /*-----------------------------------------------------------------------------
 710 Function:
 711   Find extrema between straight line C1 and hyperbola C2.
 712
 713 Method:
 714   Let P1=C1(u1) and P2=C2(u2) be two solution points
 715         D the direction of straight line C1
 716         T the tangent at point P2;
 717   Then, ( P1P2.D = 0. (1)
 718         ( P1P2.T = 0. (2)
 719   Let O1 and O2 be the origins of C1 and C2;
 720   Then, (1) <=> (O1P2-u1*D).D = 0.         as O1P1 = u1*D
 721              <=> u1 = O1P2.D               as D.D = 1.
 722          (2) <=> (P1O2 + O2P2).T= 0.
 723              <=> ((P2O1.D)D+O1O2 + O2P2).T = 0.  as P1O1 = -u1*D = (P2O1.D)D
 724              <=> (((P2O2+O2O1).D)D+O1O2 + O2P2).T = 0.
 725              <=> (P2O2.D)(D.T)+((O2O1.D)D-O2O1).T + O2P2.T= 0.
 726   We are in the reference of the hyperbola; let:
 727          by writing P (R* Chu, r* Shu, 0.0)
 728          and Chu = (v**2 + 1)/(2*v) ,
 729              Shu = (V**2 - 1)/(2*v)
 730
 731          T(R*Shu, r*Chu)
 732          D (Dx,Dy,Dz),
 733          V (Vx,Vy,Vz) = (O2O1.D)D-O2O1;
 734
 735   Then we obtain the following equation by v:
 736          (-2*R*r*Dx*Dy - R*R*Dx*Dx-r*r*Dy*Dy + R*R + r*r)     * v**4  +
 737          (2*R*Vx + 2*r*Vy)                                    * v**3  +
 738          (-2*R*Vx + 2*r*Vy)                                   * v     +
 739          (-2*R*r*Dx*Dy - (R*R*Dx*Dx-r*r*Dy*Dy + R*R + r*r))  = 0
 740
 741
 742   Use the algorithm math_DirectPolynomialRoots to solve this equation.
 743 -----------------------------------------------------------------------------*/
 744   myIsPar = Standard_False;
 745   myDone = Standard_False;
 746   myNbExt = 0;
 747
 748 // Calculate T1 in the reference of the hyperbola...
 749   gp_Dir D = C1.Direction();
 750   gp_Dir D1 = D;
 751   gp_Dir x2, y2, z2;
 752   x2 = C2.XAxis().Direction();
 753   y2 = C2.YAxis().Direction();
 754   z2 = C2.Axis().Direction();
 755   Standard_Real Dx = D.Dot(x2);
 756   Standard_Real Dy = D.Dot(y2);
 757   Standard_Real Dz = D.Dot(z2);
 758   D.SetCoord(Dx,Dy,Dz);
 759
 760 // Calculate V ...
 761   gp_Pnt O1 = C1.Location();
 762   gp_Pnt O2 = C2.Location();
 763   gp_Vec O2O1 (O2,O1);
 764   O2O1.SetCoord(O2O1.Dot(x2), O2O1.Dot(y2), O2O1.Dot(z2));
 765   gp_XYZ Vxyz = (D.XYZ()*(O2O1.Dot(D)))-O2O1.XYZ();
 766   Standard_Real Vx = Vxyz.X();
 767   Standard_Real Vy = Vxyz.Y();
 768
 769 // Calculate coefficients of the equation by v
 770   Standard_Real R = C2.MajorRadius();
 771   Standard_Real r = C2.MinorRadius();
 772   Standard_Real a = -2*R*r*Dx*Dy;
 773   Standard_Real b = -R*R*Dx*Dx - r*r*Dy*Dy + R*R + r*r;
 774   Standard_Real A1 = a + b;
 775   Standard_Real A2 = 2*R*Vx + 2*r*Vy;
 776   Standard_Real A4 = -2*R*Vx + 2*r*Vy;
 777   Standard_Real A5 = a - b;
 778
 779   math_DirectPolynomialRoots Sol(A1,A2,0.0,A4, A5);
 780   if (!Sol.IsDone()) { return; }
 781
 782 // Store solutions ...
 783   gp_Pnt P1,P2;
 784   Standard_Real U1,U2, v;
 785   Standard_Integer NbSol = Sol.NbSolutions();
 786   for (Standard_Integer NoSol = 1; NoSol <= NbSol; NoSol++) {
 787     v = Sol.Value(NoSol);
 788     if (v > 0.0) {
 789       U2 = Log(v);
 790       P2 = ElCLib::Value(U2,C2);
 791       U1 = (gp_Vec(O1,P2)).Dot(D1);
 792       P1 = ElCLib::Value(U1,C1);
 793       mySqDist[myNbExt] = P1.SquareDistance(P2);
 794       myPoint[myNbExt][0] = Extrema_POnCurv(U1,P1);
 795       myPoint[myNbExt][1] = Extrema_POnCurv(U2,P2);
 796       myNbExt++;
 797     }
 798   }
 799   myDone = Standard_True;
 800 }
 801 //=======================================================================
 802 //function : Extrema_ExtElC
 803 //purpose  :
 804 //=======================================================================
 805 Extrema_ExtElC::Extrema_ExtElC (const gp_Lin& C1,
 806                                 const gp_Parab& C2)
 807 {
 808 /*-----------------------------------------------------------------------------
 809 Function:
 810   Find extreme distances between straight line C1 and parabole C2.
 811
 812 Method:
 813    Let P1=C1(u1) and P2=C2(u2) be two solution points
 814         D the direction of straight line C1
 815         T the tangent to point P2;
 816   Then, ( P1P2.D = 0. (1)
 817         ( P1P2.T = 0. (2)
 818   Let O1 and O2 be the origins of C1 and C2;
 819   Then, (1) <=> (O1P2-u1*D).D = 0.         as O1P1 = u1*D
 820              <=> u1 = O1P2.D               as D.D = 1.
 821          (2) <=> (P1O2 + O2P2).T= 0.
 822              <=> ((P2O1.D)D+O1O2 + O2P2).T = 0.  as P1O1 = -u1*D = (P2O1.D)D
 823              <=> (((P2O2+O2O1).D)D+O1O2 + O2P2).T = 0.
 824              <=> (P2O2.D)(D.T)+((O2O1.D)D-O2O1).T + O2P2.T = 0.
 825   We are in the reference of the parabola; let:
 826          P2 (y*y/(2*p), y, 0)
 827          T (y/p, 1, 0)
 828          D (Dx,Dy,Dz),
 829          V (Vx,Vy,Vz) = (O2O1.D)D-O2O1;
 830
 831   Then, get the following equation by y:
 832      ((1-Dx*Dx)/(2*p*p))            *  y*y*y  +        A1
 833      (-3*Dx*Dy/(2*p))               *  y*y    +        A2
 834      (1-Dy*Dy + Vx/p)               *  y      +        A3
 835         Vy                          = 0.               A4
 836
 837   Use the algorithm math_DirectPolynomialRoots to solve this equation.
 838 -----------------------------------------------------------------------------*/
 839   myIsPar = Standard_False;
 840   myDone = Standard_False;
 841   myNbExt = 0;
 842
 843 // Calculate T1 in the reference of the parabola...
 844   gp_Dir D = C1.Direction();
 845   gp_Dir D1 = D;
 846   gp_Dir x2, y2, z2;
 847   x2 = C2.XAxis().Direction();
 848   y2 = C2.YAxis().Direction();
 849   z2 = C2.Axis().Direction();
 850   Standard_Real Dx = D.Dot(x2);
 851   Standard_Real Dy = D.Dot(y2);
 852   Standard_Real Dz = D.Dot(z2);
 853   D.SetCoord(Dx,Dy,Dz);
 854
 855 // Calculate V ...
 856   gp_Pnt O1 = C1.Location();
 857   gp_Pnt O2 = C2.Location();
 858   gp_Vec O2O1 (O2,O1);
 859   O2O1.SetCoord(O2O1.Dot(x2), O2O1.Dot(y2), O2O1.Dot(z2));
 860   gp_XYZ Vxyz = (D.XYZ()*(O2O1.Dot(D)))-O2O1.XYZ();
 861
 862 // Calculate coefficients of the equation by y
 863   Standard_Real P = C2.Parameter();
 864   Standard_Real A1 = (1-Dx*Dx)/(2.0*P*P);
 865   Standard_Real A2 = (-3.0*Dx*Dy/(2.0*P));
 866   Standard_Real A3 = (1 - Dy*Dy + Vxyz.X()/P);
 867   Standard_Real A4 = Vxyz.Y();
 868
 869   math_DirectPolynomialRoots Sol(A1,A2,A3,A4);
 870   if (!Sol.IsDone()) { return; }
 871
 872 // Storage of solutions ...
 873   gp_Pnt P1,P2;
 874   Standard_Real U1,U2;
 875   Standard_Integer NbSol = Sol.NbSolutions();
 876   for (Standard_Integer NoSol = 1; NoSol <= NbSol; NoSol++) {
 877     U2 = Sol.Value(NoSol);
 878     P2 = ElCLib::Value(U2,C2);
 879     U1 = (gp_Vec(O1,P2)).Dot(D1);
 880     P1 = ElCLib::Value(U1,C1);
 881     mySqDist[myNbExt] = P1.SquareDistance(P2);
 882     myPoint[myNbExt][0] = Extrema_POnCurv(U1,P1);
 883     myPoint[myNbExt][1] = Extrema_POnCurv(U2,P2);
 884     myNbExt++;
 885   }
 886   myDone = Standard_True;
 887 }
 888 //=======================================================================
 889 //function : Extrema_ExtElC
 890 //purpose  :
 891 //=======================================================================
 892 Extrema_ExtElC::Extrema_ExtElC (const gp_Circ& C1,
 893                                 const gp_Circ& C2)
 894 {
 895   Standard_Boolean bIsSamePlane, bIsSameAxe;
 896   Standard_Real aTolD, aTolD2, aTolA, aD2, aDC2;
 897   gp_Pnt aPc1, aPc2;
 898   gp_Dir aDc1, aDc2;
 899   //
 900   myIsPar = Standard_False;
 901   myDone = Standard_False;
 902   myNbExt = 0;
 903   //
 904   aTolA=Precision::Angular();
 905   aTolD=Precision::Confusion();
 906   aTolD2=aTolD*aTolD;
 907   //
 908   aPc1=C1.Location();
 909   aDc1=C1.Axis().Direction();
 910
 911   aPc2=C2.Location();
 912   aDc2=C2.Axis().Direction();
 913   gp_Pln aPlc1(aPc1, aDc1);
 914   //
 915   aD2=aPlc1.SquareDistance(aPc2);
 916   bIsSamePlane=aDc1.IsParallel(aDc2, aTolA) && aD2<aTolD2;
 917   if (!bIsSamePlane) {
 918     return;
 919   }
 920
 921   // Here, both circles are in the same plane.
 922
 923   //
 924   aDC2=aPc1.SquareDistance(aPc2);
 925   bIsSameAxe=aDC2<aTolD2;
 926   //
 927   if(bIsSameAxe)
 928   {
 929     myIsPar = Standard_True;
 930     myNbExt = 1;
 931     myDone = Standard_True;
 932     const Standard_Real aDR = C1.Radius() - C2.Radius();
 933     mySqDist[0] = aDR*aDR;
 934     return;
 935   }
 936
 937   Standard_Boolean bIn, bOut;
 938   Standard_Integer j1, j2;
 939   Standard_Real aR1, aR2, aD12, aT11, aT12, aT21, aT22;
 940   gp_Circ aC1, aC2;
 941   gp_Pnt aP11, aP12, aP21, aP22;
 942   //
 943   myDone = Standard_True;
 944   //
 945   aR1 = C1.Radius();
 946   aR2 = C2.Radius();
 947   //
 948   j1 = 0;
 949   j2 = 1;
 950   aC1 = C1;
 951   aC2 = C2;
 952   if (aR2 > aR1)
 953   {
 954     j1 = 1;
 955     j2 = 0;
 956     aC1 = C2;
 957     aC2 = C1;
 958   }
 959   //
 960   aR1 = aC1.Radius(); // max radius
 961   aR2 = aC2.Radius(); // min radius
 962   //
 963   aPc1 = aC1.Location();
 964   aPc2 = aC2.Location();
 965   //
 966   aD12 = aPc1.Distance(aPc2);
 967   gp_Vec aVec12(aPc1, aPc2);
 968   gp_Dir aDir12(aVec12);
 969   //
 970   // 1. Four common solutions
 971   myNbExt = 4;
 972   //
 973   aP11.SetXYZ(aPc1.XYZ() - aR1*aDir12.XYZ());
 974   aP12.SetXYZ(aPc1.XYZ() + aR1*aDir12.XYZ());
 975   aP21.SetXYZ(aPc2.XYZ() - aR2*aDir12.XYZ());
 976   aP22.SetXYZ(aPc2.XYZ() + aR2*aDir12.XYZ());
 977   //
 978   aT11 = ElCLib::Parameter(aC1, aP11);
 979   aT12 = ElCLib::Parameter(aC1, aP12);
 980   aT21 = ElCLib::Parameter(aC2, aP21);
 981   aT22 = ElCLib::Parameter(aC2, aP22);
 982   //
 983   // P11, P21
 984   myPoint[0][j1].SetValues(aT11, aP11);
 985   myPoint[0][j2].SetValues(aT21, aP21);
 986   mySqDist[0] = aP11.SquareDistance(aP21);
 987   // P11, P22
 988   myPoint[1][j1].SetValues(aT11, aP11);
 989   myPoint[1][j2].SetValues(aT22, aP22);
 990   mySqDist[1] = aP11.SquareDistance(aP22);
 991   //
 992   // P12, P21
 993   myPoint[2][j1].SetValues(aT12, aP12);
 994   myPoint[2][j2].SetValues(aT21, aP21);
 995   mySqDist[2] = aP12.SquareDistance(aP21);
 996   //
 997   // P12, P22
 998   myPoint[3][j1].SetValues(aT12, aP12);
 999   myPoint[3][j2].SetValues(aT22, aP22);
1000   mySqDist[3] = aP12.SquareDistance(aP22);
1001   //
1002   // 2. Check for intersections
1003   bOut = aD12 > (aR1 + aR2 + aTolD);
1004   bIn = aD12 < (aR1 - aR2 - aTolD);
1005   if (!bOut && !bIn)
1006   {
1007     Standard_Boolean bNbExt6;
1008     Standard_Real aAlpha, aBeta, aT[2], aVal, aDist2;
1009     gp_Pnt aPt, aPL1, aPL2;
1010     gp_Dir aDLt;
1011     //
1012     aAlpha = 0.5*(aR1*aR1 - aR2*aR2 + aD12*aD12) / aD12;
1013     aVal = aR1*aR1 - aAlpha*aAlpha;
1014     if (aVal < 0.)
1015     {// see pkv/900/L4 for details
1016       aVal = -aVal;
1017     }
1018     aBeta = Sqrt(aVal);
1019     //aBeta=Sqrt(aR1*aR1-aAlpha*aAlpha);
1020     //--
1021     aPt.SetXYZ(aPc1.XYZ() + aAlpha*aDir12.XYZ());
1022     //
1023     aDLt = aDc1^aDir12;
1024     aPL1.SetXYZ(aPt.XYZ() + aBeta*aDLt.XYZ());
1025     aPL2.SetXYZ(aPt.XYZ() - aBeta*aDLt.XYZ());
1026     //
1027     aDist2 = aPL1.SquareDistance(aPL2);
1028     bNbExt6 = aDist2 > aTolD2;
1029     //
1030     myNbExt = 5;// just in case. see pkv/900/L4 for details
1031     aT[j1] = ElCLib::Parameter(aC1, aPL1);
1032     aT[j2] = ElCLib::Parameter(aC2, aPL1);
1033     myPoint[4][j1].SetValues(aT[j1], aPL1);
1034     myPoint[4][j2].SetValues(aT[j2], aPL1);
1035     mySqDist[4] = 0.;
1036     //
1037     if (bNbExt6)
1038     {
1039       myNbExt = 6;
1040       aT[j1] = ElCLib::Parameter(aC1, aPL2);
1041       aT[j2] = ElCLib::Parameter(aC2, aPL2);
1042       myPoint[5][j1].SetValues(aT[j1], aPL2);
1043       myPoint[5][j2].SetValues(aT[j2], aPL2);
1044       mySqDist[5] = 0.;
1045     }
1046     //
1047   }// if (!bOut || !bIn) {
1048 }
1049
1050 //=======================================================================
1051 //function : IsDone
1052 //purpose  :
1053 //=======================================================================
1054 Standard_Boolean Extrema_ExtElC::IsDone () const {
1055   return myDone;
1056 }
1057 //=======================================================================
1058 //function : IsParallel
1059 //purpose  :
1060 //=======================================================================
1061 Standard_Boolean Extrema_ExtElC::IsParallel () const
1062 {
1063   if (!IsDone()) {
1064     throw StdFail_NotDone();
1065   }
1066   return myIsPar;
1067 }
1068 //=======================================================================
1069 //function : NbExt
1070 //purpose  :
1071 //=======================================================================
1072 Standard_Integer Extrema_ExtElC::NbExt () const
1073 {
1074   if (!IsDone())
1075   {
1076     throw StdFail_NotDone();
1077   }
1078   return myNbExt;
1079 }
1080 //=======================================================================
1081 //function : SquareDistance
1082 //purpose  :
1083 //=======================================================================
1084 Standard_Real Extrema_ExtElC::SquareDistance (const Standard_Integer N) const
1085 {
1086   if (N < 1 || N > NbExt())
1087   {
1088     throw Standard_OutOfRange();
1089   }
1090
1091   return mySqDist[N - 1];
1092 }
1093 //=======================================================================
1094 //function : Points
1095 //purpose  :
1096 //=======================================================================
1097 void Extrema_ExtElC::Points (const Standard_Integer N,
1098                              Extrema_POnCurv& P1,
1099                              Extrema_POnCurv& P2) const
1100 {
1101   if (N < 1 || N > NbExt())
1102   {
1103     throw Standard_OutOfRange();
1104   }
1105
1106   P1 = myPoint[N-1][0];
1107   P2 = myPoint[N-1][1];
1108 }
1109
1110
1111 //=======================================================================
1112 //function : RefineDir
1113 //purpose  :
1114 //=======================================================================
1115 void RefineDir(gp_Dir& aDir)
1116 {
1117   Standard_Integer i, j, k, iK;
1118   Standard_Real aCx[3], aEps, aX1, aX2, aOne;
1119   //
1120   iK=3;
1121   aEps=RealEpsilon();
1122   aDir.Coord(aCx[0], aCx[1], aCx[2]);
1123   //
1124   for (i=0; i<iK; ++i) {
1125     aOne=(aCx[i]>0.) ? 1. : -1.;
1126     aX1=aOne-aEps;
1127     aX2=aOne+aEps;
1128     //
1129     if (aCx[i]>aX1 && aCx[i]<aX2) {
1130       j=(i+1)%iK;
1131       k=(i+2)%iK;
1132       aCx[i]=aOne;
1133       aCx[j]=0.;
1134       aCx[k]=0.;
1135       aDir.SetCoord(aCx[0], aCx[1], aCx[2]);
1136       return;
1137     }
1138   }
1139 }