src/Extrema/Extrema_ExtElC.cxx

   1
   2 #include <stdio.h>
   3
   4 #include <Extrema_ExtElC.ixx>
   5 #include <StdFail_InfiniteSolutions.hxx>
   6 #include <StdFail_NotDone.hxx>
   7 #include <ElCLib.hxx>
   8 #include <math_TrigonometricFunctionRoots.hxx>
   9 #include <math_DirectPolynomialRoots.hxx>
  10 #include <Standard_OutOfRange.hxx>
  11 #include <Standard_NotImplemented.hxx>
  12 #include <Precision.hxx>
  13 #include <Extrema_ExtPElC.hxx>
  14 #include <IntAna_QuadQuadGeo.hxx>
  15 #include <Extrema_ExtPElC.hxx>
  16
  17 #include <gp_Pln.hxx>
  18 #include <gp_Ax3.hxx>
  19 #include <gp_Ax2.hxx>
  20 #include <gp_Pnt.hxx>
  21 #include <gp_Dir.hxx>
  22 #include <gp_Ax1.hxx>
  23
  24
  25
  26 //======================================================================
  27 //==  Classe Interne (Donne des racines classees d un polynome trigo)
  28 //==  Code duplique avec IntAna_IntQuadQuad.cxx (lbr le 26 mars 98)
  29 //==  Solution fiable aux problemes de coefficients proches de 0
  30 //==  avec essai de rattrapage si coeff<1.e-10 (jct le 27 avril 98)
  31 //======================================================================
  32 static const Standard_Real PIpPI=Standard_PI+Standard_PI;
  33
  34 class ExtremaExtElC_TrigonometricRoots {
  35  private:
  36   Standard_Real Roots[4];
  37   Standard_Boolean done;
  38   Standard_Integer NbRoots;
  39   Standard_Boolean infinite_roots;
  40  public:
  41   ExtremaExtElC_TrigonometricRoots(const Standard_Real CC
  42                      ,const Standard_Real SC
  43                      ,const Standard_Real C
  44                      ,const Standard_Real S
  45                      ,const Standard_Real Cte
  46                      ,const Standard_Real Binf
  47                      ,const Standard_Real Bsup);
  48
  49   Standard_Boolean IsDone() { return(done); }
  50
  51   Standard_Boolean IsARoot(Standard_Real u) {
  52     for(Standard_Integer i=0 ; i<NbRoots; i++) {
  53       if(Abs(u - Roots[i])<=RealEpsilon()) return(Standard_True);
  54       if(Abs(u - Roots[i]-PIpPI)<=RealEpsilon()) return(Standard_True);
  55     }
  56     return(Standard_False);
  57   }
  58
  59   Standard_Integer NbSolutions() {
  60     if(!done) StdFail_NotDone::Raise();
  61     return(NbRoots);
  62   }
  63   Standard_Boolean InfiniteRoots() {
  64     if(!done) StdFail_NotDone::Raise();
  65     return(infinite_roots);
  66   }
  67   Standard_Real    Value(const Standard_Integer& n) {
  68     if((!done)||(n>NbRoots)) StdFail_NotDone::Raise();
  69     return(Roots[n-1]);
  70   }
  71 };
  72 //----------------------------------------------------------------------
  73 ExtremaExtElC_TrigonometricRoots::ExtremaExtElC_TrigonometricRoots(const Standard_Real CC
  74                                        ,const Standard_Real SC
  75                                        ,const Standard_Real C
  76                                        ,const Standard_Real S
  77                                        ,const Standard_Real Cte
  78                                        ,const Standard_Real Binf
  79                                        ,const Standard_Real Bsup) {
  80
  81   Standard_Integer i ;
  82   Standard_Integer nbessai = 1;
  83   Standard_Real cc = CC, sc = SC, c = C, s = S, cte = Cte;
  84   done=Standard_False;
  85   while (nbessai<=2 && !done) {
  86     //-- F= AA*CN*CN+2*BB*CN*SN+CC*CN+DD*SN+EE;
  87     math_TrigonometricFunctionRoots MTFR(cc,sc,c,s,cte,Binf,Bsup);
  88     if(MTFR.IsDone()) {
  89       done=Standard_True;
  90       if(MTFR.InfiniteRoots()) {
  91         infinite_roots=Standard_True;
  92       }
  93       else {
  94         NbRoots=MTFR.NbSolutions();
  95         for( i=0;i<NbRoots;i++) {
  96           Roots[i]=MTFR.Value(i+1);
  97           if(Roots[i]<0.0) Roots[i]+=PI+PI;
  98           if(Roots[i]>(PI+PI)) Roots[i]-=PI+PI;
  99         }
 100         Standard_Boolean Triee;
 101       Standard_Integer j;
 102
 103         //-- La recherche directe donne n importe quoi.
 104         // modified by OCC  Tue Oct  3 18:41:27 2006.BEGIN
 105         Standard_Real aMaxCoef = Max(CC,SC);
 106         aMaxCoef = Max(aMaxCoef,C);
 107         aMaxCoef = Max(aMaxCoef,S);
 108         aMaxCoef = Max(aMaxCoef,Cte);
 109         const Standard_Real aPrecision = Max(1.e-8,1.e-12*aMaxCoef);
 110         // modified by OCC  Tue Oct  3 18:41:33 2006.END
 111
 112         Standard_Integer SvNbRoots=NbRoots;
 113         for(i=0;i<SvNbRoots;i++) {
 114           Standard_Real y;
 115         Standard_Real co=cos(Roots[i]);
 116           Standard_Real si=sin(Roots[i]);
 117           y=co*(CC*co + (SC+SC)*si + C) + S*si + Cte;
 118           // modified by OCC  Tue Oct  3 18:43:00 2006
 119           if(Abs(y)>aPrecision) {
 120 #ifdef DEB
 121             printf("\n**IntAna_IntQuadQuad** Solution : %g ( %g cos2 + 2  %g cos sin + %g cos + %g sin  + %g) = %g\n",
 122                    Roots[i],CC,SC,C,S,Cte,y);
 123 #endif
 124             NbRoots--;
 125             Roots[i]=1000.0;
 126           }
 127         }
 128
 129         do {
 130           Triee=Standard_True;
 131           for(i=1,j=0;i<SvNbRoots;i++,j++) {
 132             if(Roots[i]<Roots[j]) {
 133               Triee=Standard_False;
 134               Standard_Real t=Roots[i]; Roots[i]=Roots[j]; Roots[j]=t;
 135             }
 136           }
 137         }
 138         while(!Triee);
 139         infinite_roots=Standard_False;
 140         if(NbRoots==0) { //--!!!!! Detection du cas Pol = Cte ( 1e-50 ) !!!!
 141           if((Abs(CC) + Abs(SC) + Abs(C) + Abs(S)) < 1e-10) {
 142             if(Abs(Cte) < 1e-10)  {
 143               infinite_roots=Standard_True;
 144             }
 145           }
 146         }
 147       }
 148     }
 149     else {
 150         // on essaie en mettant les tres petits coeff. a ZERO
 151       if (Abs(CC)<1e-10) cc = 0.0;
 152       if (Abs(SC)<1e-10) sc = 0.0;
 153       if (Abs(C)<1e-10) c = 0.0;
 154       if (Abs(S)<1e-10) s = 0.0;
 155       if (Abs(Cte)<1e-10) cte = 0.0;
 156       nbessai++;
 157     }
 158   }
 159 }
 160
 161 //=============================================================================
 162
 163 Extrema_ExtElC::Extrema_ExtElC () { myDone = Standard_False; }
 164 //=============================================================================
 165
 166 Extrema_ExtElC::Extrema_ExtElC (const gp_Lin& C1, const gp_Lin& C2,
 167                                 const Standard_Real)
 168 // Fonction:
 169 //    Recherche de la distance minimale entre 2 droites.
 170
 171 // Methode:
 172 //   Soit D1 et D2, les 2 directions des droites C1 et C2.
 173 //   2 cas sont consideres:
 174 //   1- si Angle(D1,D2) < AngTol, les droites sont paralleles.
 175 //      La distance est la distance entre un point quelconque de C1 et la droite
 176 //      C2.
 177 //   2- si Angle(D1,D2) > AngTol:
 178 //      Soit P1=C1(u1) et P2=C2(u2) les 2 points solutions:
 179 //      Alors, ( P1P2.D1 = 0. (1)
 180 //             ( P1P2.D2 = 0. (2)
 181 //   Soit O1 et O2 les origines de C1 et C2;
 182 //   Alors, (1) <=> (O1P2-u1*D1).D1 = 0.       car O1P1 = u1*D1
 183 //           <=> u1 = O1P2.D1               car D1.D1 = 1.
 184 //          (2) <=> (P1O2+u2*D2).D2 = 0.       car O2P2 = u2*D2
 185 //           <=> u2 = O2P1.D2               car D2.D2 = 1.
 186 //           <=> u2 = (O2O1+O1P1).D2
 187 //           <=> u2 = O2O1.D2+((O1P2.T1)T1).T2) car O1P1 = u1*T1 = (O1P2.T1)T1
 188 //           <=> u2 = O2O1.D2+(((O1O2+O2P2).D1)D1).D2)
 189 //           <=> u2 = O2O1.D2+((O1O2.D1)D1).D2)+(O2P2.D1)(D1.D2)
 190 //           <=> u2 = ((O1O2.D1)D1-O1O2).D2 + u2*(D2.D1)(D1.D2)
 191 //           <=> u2 = (((O1O2.D1)D1-O1O2).D2) / 1.-(D1.D2)**2
 192 {
 193   myDone = Standard_False;
 194   myNbExt = 0;
 195
 196   gp_Dir D1 = C1.Position().Direction();
 197   gp_Dir D2 = C2.Position().Direction();
 198   // MSV 16/01/2000: avoid "divide by zero"
 199   Standard_Real D1DotD2 = D1.Dot(D2);
 200   Standard_Real aSin = 1.-D1DotD2*D1DotD2;
 201   if (aSin < gp::Resolution() ||
 202       D1.IsParallel(D2, Precision::Angular())) {
 203     myIsPar = Standard_True;
 204     mySqDist[0] = C2.SquareDistance(C1.Location());
 205     mySqDist[1] = mySqDist[0];
 206   }
 207   else {
 208     myIsPar = Standard_False;
 209     gp_Pnt O1 = C1.Location();
 210     gp_Pnt O2 = C2.Location();
 211     gp_Vec O1O2 (O1,O2);
 212     Standard_Real U2 = (D1.XYZ()*(O1O2.Dot(D1))-(O1O2.XYZ())).Dot(D2.XYZ());
 213     if( Precision::IsInfinite(U2) ) {
 214       myIsPar = Standard_True;
 215       mySqDist[0] = C2.SquareDistance(C1.Location());
 216       mySqDist[1] = mySqDist[0];
 217     }
 218     else {
 219       U2 /= aSin;
 220       if( Precision::IsInfinite(U2) ) {
 221         myIsPar = Standard_True;
 222         mySqDist[0] = C2.SquareDistance(C1.Location());
 223     mySqDist[1] = mySqDist[0];
 224       }
 225       else {
 226         gp_Pnt P2(ElCLib::Value(U2,C2));
 227         Standard_Real U1 = (gp_Vec(O1,P2)).Dot(D1);
 228         if( Precision::IsInfinite(U1) ) {
 229           myIsPar = Standard_True;
 230           mySqDist[0] = C2.SquareDistance(C1.Location());
 231       mySqDist[1] = mySqDist[0];
 232         }
 233         else {
 234           gp_Pnt P1(ElCLib::Value(U1,C1));
 235           mySqDist[myNbExt] = P1.SquareDistance(P2);
 236           myPoint[myNbExt][0] = Extrema_POnCurv(U1,P1);
 237           myPoint[myNbExt][1] = Extrema_POnCurv(U2,P2);
 238           myNbExt = 1;
 239         }
 240       }
 241     }
 242   }
 243   myDone = Standard_True;
 244 }
 245 //=============================================================================
 246
 247 Extrema_ExtElC::Extrema_ExtElC (const gp_Lin& C1, const gp_Circ& C2,
 248                                 const Standard_Real)
 249 /*-----------------------------------------------------------------------------
 250 Fonction:
 251    Recherche des distances extremales entre la droite C1 et le cercle C2.
 252
 253 Methode:
 254    Soit P1=C1(u1) et P2=C2(u2) deux points solutions
 255         D la direction de la droite C1
 256         T la tangente au point P2;
 257   Alors, ( P1P2.D = 0. (1)
 258          ( P1P2.T = 0. (2)
 259   Soit O1 et O2 les origines de C1 et C2;
 260   Alors, (1) <=> (O1P2-u1*D).D = 0.         car O1P1 = u1*D
 261              <=> u1 = O1P2.D                car D.D = 1.
 262          (2) <=> P1O2.T = 0.                car O2P2.T = 0.
 263              <=> ((P2O1.D)D+O1O2).T = 0.    car P1O1 = -u1*D = (P2O1.D)D
 264              <=> (((P2O2+O2O1).D)D+O1O2).T = 0.
 265              <=> ((P2O2.D)(D.T)+((O2O1.D)D-O2O1).T = 0.
 266   On se place dans le repere du cercle; soit:
 267          Cos = Cos(u2) et Sin = Sin(u2),
 268          P2 (R*Cos,R*Sin,0.),
 269          T (-R*Sin,R*Cos,0.),
 270          D (Dx,Dy,Dz),
 271          V (Vx,Vy,Vz) = (O2O1.D)D-O2O1;
 272   Alors, on obtient l'equation en Cos et Sin suivante:
 273     -(2*R*R*Dx*Dy)   * Cos**2  +       A1
 274    R*R*(Dx**2-Dy**2) * Cos*Sin +    2* A2
 275          R*Vy        * Cos     +       A3
 276         -R*Vx        * Sin     +       A4
 277       R*R*Dx*Dy                = 0.    A5
 278      On utilise l'algorithme math_TrigonometricFunctionRoots pour resoudre
 279     cette equation.
 280 -----------------------------------------------------------------------------*/
 281 {
 282   myIsPar = Standard_False;
 283   myDone = Standard_False;
 284   myNbExt = 0;
 285
 286 // Calcul de T1 dans le repere du cercle ...
 287   gp_Dir D = C1.Direction();
 288   gp_Dir D1 = D;
 289   gp_Dir x2, y2, z2;
 290   x2 = C2.XAxis().Direction();
 291   y2 = C2.YAxis().Direction();
 292   z2 = C2.Axis().Direction();
 293   Standard_Real Dx = D.Dot(x2);
 294   Standard_Real Dy = D.Dot(y2);
 295   Standard_Real Dz = D.Dot(z2);
 296   D.SetCoord(Dx,Dy,Dz);
 297
 298 // Calcul de V dans le repere du cercle:
 299   gp_Pnt O1 = C1.Location();
 300   gp_Pnt O2 = C2.Location();
 301   gp_Vec O2O1 (O2,O1);
 302   O2O1.SetCoord(O2O1.Dot(x2), O2O1.Dot(y2), O2O1.Dot(z2));
 303   gp_XYZ Vxyz = (D.XYZ()*(O2O1.Dot(D)))-O2O1.XYZ();
 304
 305 // Calcul des coefficients de l equation en Cos et Sin ...
 306   Standard_Real R = C2.Radius();
 307   Standard_Real A5 = R*R*Dx*Dy;
 308   Standard_Real A1 = -2.*A5;
 309   Standard_Real A2 = R*R*(Dx*Dx-Dy*Dy)/2.0;
 310   Standard_Real A3 = R*Vxyz.Y();
 311   Standard_Real A4 = -R*Vxyz.X();
 312  // Standard_Real TolU = Tol * R;
 313
 314
 315   if(fabs(A5) <= 1.e-12) A5 = 0.;
 316   if(fabs(A1) <= 1.e-12) A1 = 0.;
 317   if(fabs(A2) <= 1.e-12) A2 = 0.;
 318   if(fabs(A3) <= 1.e-12) A3 = 0.;
 319   if(fabs(A4) <= 1.e-12) A4 = 0.;
 320
 321
 322   ExtremaExtElC_TrigonometricRoots Sol(A1,A2,A3,A4,A5,0.,PI+PI);
 323   if (!Sol.IsDone()) { return; }
 324   if (Sol.InfiniteRoots()) {
 325     myIsPar = Standard_True;
 326     mySqDist[0] = R*R;
 327     myDone = Standard_True;
 328     return;
 329   }
 330 // Stockage des solutions ...
 331   gp_Pnt P1,P2;
 332   Standard_Real U1,U2;
 333   Standard_Integer NbSol = Sol.NbSolutions();
 334   for (Standard_Integer NoSol = 1; NoSol <= NbSol; NoSol++) {
 335     U2 = Sol.Value(NoSol);
 336     P2 = ElCLib::Value(U2,C2);
 337     U1 = (gp_Vec(O1,P2)).Dot(D1);
 338     P1 = ElCLib::Value(U1,C1);
 339     mySqDist[myNbExt] = P1.SquareDistance(P2);
 340     myPoint[myNbExt][0] = Extrema_POnCurv(U1,P1);
 341     myPoint[myNbExt][1] = Extrema_POnCurv(U2,P2);
 342     myNbExt++;
 343   }
 344   myDone = Standard_True;
 345 }
 346
 347
 348
 349 Extrema_ExtElC::Extrema_ExtElC (const gp_Lin& C1, const gp_Elips& C2)
 350 {
 351 /*-----------------------------------------------------------------------------
 352 Fonction:
 353    Recherche des distances extremales entre la droite C1 et l ellipse C2.
 354
 355 Methode:
 356    Soit P1=C1(u1) et P2=C2(u2) deux points solutions
 357         D la direction de la droite C1
 358         T la tangente au point P2;
 359   Alors, ( P1P2.D = 0. (1)
 360          ( P1P2.T = 0. (2)
 361   Soit O1 et O2 les origines de C1 et C2;
 362   Alors, (1) <=> (O1P2-u1*D).D = 0.         car O1P1 = u1*D
 363              <=> u1 = O1P2.D                car D.D = 1.
 364          (2) <=> P1O2.T = 0.                car O2P2.T = 0.
 365              <=> ((P2O1.D)D+O1O2).T = 0.    car P1O1 = -u1*D = (P2O1.D)D
 366              <=> (((P2O2+O2O1).D)D+O1O2).T = 0.
 367              <=> ((P2O2.D)(D.T)+((O2O1.D)D-O2O1).T = 0.
 368   On se place dans le repere de l ellipse; soit:
 369          Cos = Cos(u2) et Sin = Sin(u2),
 370          P2 (MajR*Cos,MinR*Sin,0.),
 371          T (-MajR*Sin,MinR*Cos,0.),
 372          D (Dx,Dy,Dz),
 373          V (Vx,Vy,Vz) = (O2O1.D)D-O2O1;
 374   Alors, on obtient l'equation en Cos et Sin suivante:
 375     -(2*MajR*MinR*Dx*Dy)             * Cos**2  +
 376    (MajR*MajR*Dx**2-MinR*MinR*Dy**2) * Cos*Sin +
 377          MinR*Vy                     * Cos     +
 378        - MajR*Vx                     * Sin     +
 379       MinR*MajR*Dx*Dy                = 0.
 380      On utilise l'algorithme math_TrigonometricFunctionRoots pour resoudre
 381     cette equation.
 382 -----------------------------------------------------------------------------*/
 383   myIsPar = Standard_False;
 384   myDone = Standard_False;
 385   myNbExt = 0;
 386
 387 // Calcul de T1 dans le repere de l'ellipse ...
 388   gp_Dir D = C1.Direction();
 389   gp_Dir D1 = D;
 390   gp_Dir x2, y2, z2;
 391   x2 = C2.XAxis().Direction();
 392   y2 = C2.YAxis().Direction();
 393   z2 = C2.Axis().Direction();
 394   Standard_Real Dx = D.Dot(x2);
 395   Standard_Real Dy = D.Dot(y2);
 396   Standard_Real Dz = D.Dot(z2);
 397   D.SetCoord(Dx,Dy,Dz);
 398
 399 // Calcul de V ...
 400   gp_Pnt O1 = C1.Location();
 401   gp_Pnt O2 = C2.Location();
 402   gp_Vec O2O1 (O2,O1);
 403   O2O1.SetCoord(O2O1.Dot(x2), O2O1.Dot(y2), O2O1.Dot(z2));
 404   gp_XYZ Vxyz = (D.XYZ()*(O2O1.Dot(D)))-O2O1.XYZ();
 405
 406 // Calcul des coefficients de l equation en Cos et Sin ...
 407   Standard_Real MajR = C2.MajorRadius();
 408   Standard_Real MinR = C2.MinorRadius();
 409   Standard_Real A5 = MajR*MinR*Dx*Dy;
 410   Standard_Real A1 = -2.*A5;
 411   Standard_Real R2 = MajR*MajR;
 412   Standard_Real r2 = MinR*MinR;
 413   Standard_Real A2 =(R2*Dx*Dx -r2*Dy*Dy -R2 +r2)/2.0;
 414   Standard_Real A3 = MinR*Vxyz.Y();
 415   Standard_Real A4 = -MajR*Vxyz.X();
 416
 417   ExtremaExtElC_TrigonometricRoots Sol(A1,A2,A3,A4,A5,0.,PI+PI);
 418   if (!Sol.IsDone()) { return; }
 419
 420 // Stockage des solutions ...
 421   gp_Pnt P1,P2;
 422   Standard_Real U1,U2;
 423   Standard_Integer NbSol = Sol.NbSolutions();
 424   for (Standard_Integer NoSol = 1; NoSol <= NbSol; NoSol++) {
 425     U2 = Sol.Value(NoSol);
 426     P2 = ElCLib::Value(U2,C2);
 427     U1 = (gp_Vec(O1,P2)).Dot(D1);
 428     P1 = ElCLib::Value(U1,C1);
 429     mySqDist[myNbExt] = P1.SquareDistance(P2);
 430     myPoint[myNbExt][0] = Extrema_POnCurv(U1,P1);
 431     myPoint[myNbExt][1] = Extrema_POnCurv(U2,P2);
 432     myNbExt++;
 433   }
 434   myDone = Standard_True;
 435 }
 436 //=============================================================================
 437
 438 Extrema_ExtElC::Extrema_ExtElC (const gp_Lin& C1, const gp_Hypr& C2)
 439 {
 440 /*-----------------------------------------------------------------------------
 441 Fonction:
 442    Recherche des distances extremales entre la droite C1 et l'hyperbole C2.
 443
 444 Methode:
 445    Soit P1=C1(u1) et P2=C2(u2) deux points solutions
 446         D la direction de la droite C1
 447         T la tangente au point P2;
 448   Alors, ( P1P2.D = 0. (1)
 449          ( P1P2.T = 0. (2)
 450   Soit O1 et O2 les origines de C1 et C2;
 451   Alors, (1) <=> (O1P2-u1*D).D = 0.         car O1P1 = u1*D
 452              <=> u1 = O1P2.D                car D.D = 1.
 453          (2) <=> (P1O2 + O2P2).T= 0.
 454              <=> ((P2O1.D)D+O1O2 + O2P2).T = 0.  car P1O1 = -u1*D = (P2O1.D)D
 455              <=> (((P2O2+O2O1).D)D+O1O2 + O2P2).T = 0.
 456              <=> (P2O2.D)(D.T)+((O2O1.D)D-O2O1).T + O2P2.T= 0.
 457   On se place dans le repere de l'hyperbole; soit:
 458          en ecrivant P (R* Chu, r* Shu, 0.0)
 459          et Chu = (v**2 + 1)/(2*v) ,
 460             Shu = (V**2 - 1)/(2*v)
 461
 462          T(R*Shu, r*Chu)
 463          D (Dx,Dy,Dz),
 464          V (Vx,Vy,Vz) = (O2O1.D)D-O2O1;
 465
 466   Alors, on obtient l'equation en v suivante:
 467          (-2*R*r*Dx*Dy - R*R*Dx*Dx-r*r*Dy*Dy + R*R + r*r)     * v**4  +
 468          (2*R*Vx + 2*r*Vy)                                    * v**3  +
 469          (-2*R*Vx + 2*r*Vy)                                   * v     +
 470          (-2*R*r*Dx*Dy - (R*R*Dx*Dx-r*r*Dy*Dy + R*R + r*r))  = 0
 471
 472
 473      On utilise l'algorithme math_DirectPolynomialRoots pour resoudre
 474     cette equation.
 475 -----------------------------------------------------------------------------*/
 476   myIsPar = Standard_False;
 477   myDone = Standard_False;
 478   myNbExt = 0;
 479
 480 // Calcul de T1 dans le repere de l'hyperbole ...
 481   gp_Dir D = C1.Direction();
 482   gp_Dir D1 = D;
 483   gp_Dir x2, y2, z2;
 484   x2 = C2.XAxis().Direction();
 485   y2 = C2.YAxis().Direction();
 486   z2 = C2.Axis().Direction();
 487   Standard_Real Dx = D.Dot(x2);
 488   Standard_Real Dy = D.Dot(y2);
 489   Standard_Real Dz = D.Dot(z2);
 490   D.SetCoord(Dx,Dy,Dz);
 491
 492 // Calcul de V ...
 493   gp_Pnt O1 = C1.Location();
 494   gp_Pnt O2 = C2.Location();
 495   gp_Vec O2O1 (O2,O1);
 496   O2O1.SetCoord(O2O1.Dot(x2), O2O1.Dot(y2), O2O1.Dot(z2));
 497   gp_XYZ Vxyz = (D.XYZ()*(O2O1.Dot(D)))-O2O1.XYZ();
 498   Standard_Real Vx = Vxyz.X();
 499   Standard_Real Vy = Vxyz.Y();
 500
 501 // Calcul des coefficients de l equation en v
 502   Standard_Real R = C2.MajorRadius();
 503   Standard_Real r = C2.MinorRadius();
 504   Standard_Real a = -2*R*r*Dx*Dy;
 505   Standard_Real b = -R*R*Dx*Dx - r*r*Dy*Dy + R*R + r*r;
 506   Standard_Real A1 = a + b;
 507   Standard_Real A2 = 2*R*Vx + 2*r*Vy;
 508   Standard_Real A4 = -2*R*Vx + 2*r*Vy;
 509   Standard_Real A5 = a - b;
 510
 511   math_DirectPolynomialRoots Sol(A1,A2,0.0,A4, A5);
 512   if (!Sol.IsDone()) { return; }
 513
 514 // Stockage des solutions ...
 515   gp_Pnt P1,P2;
 516   Standard_Real U1,U2, v;
 517   Standard_Integer NbSol = Sol.NbSolutions();
 518   for (Standard_Integer NoSol = 1; NoSol <= NbSol; NoSol++) {
 519     v = Sol.Value(NoSol);
 520     if (v > 0.0) {
 521       U2 = Log(v);
 522       P2 = ElCLib::Value(U2,C2);
 523       U1 = (gp_Vec(O1,P2)).Dot(D1);
 524       P1 = ElCLib::Value(U1,C1);
 525       mySqDist[myNbExt] = P1.SquareDistance(P2);
 526       myPoint[myNbExt][0] = Extrema_POnCurv(U1,P1);
 527       myPoint[myNbExt][1] = Extrema_POnCurv(U2,P2);
 528       myNbExt++;
 529     }
 530   }
 531   myDone = Standard_True;
 532 }
 533 //=============================================================================
 534
 535 Extrema_ExtElC::Extrema_ExtElC (const gp_Lin& C1, const gp_Parab& C2)
 536 {
 537 /*-----------------------------------------------------------------------------
 538 Fonction:
 539    Recherche des distances extremales entre la droite C1 et la parabole C2.
 540
 541 Methode:
 542    Soit P1=C1(u1) et P2=C2(u2) deux points solutions
 543         D la direction de la droite C1
 544         T la tangente au point P2;
 545   Alors, ( P1P2.D = 0. (1)
 546          ( P1P2.T = 0. (2)
 547   Soit O1 et O2 les origines de C1 et C2;
 548   Alors, (1) <=> (O1P2-u1*D).D = 0.         car O1P1 = u1*D
 549              <=> u1 = O1P2.D                car D.D = 1.
 550          (2) <=> (P1O2 + O2P2).T= 0.
 551              <=> ((P2O1.D)D+O1O2 + O2P2).T = 0.  car P1O1 = -u1*D = (P2O1.D)D
 552              <=> (((P2O2+O2O1).D)D+O1O2 + O2P2).T = 0.
 553              <=> (P2O2.D)(D.T)+((O2O1.D)D-O2O1).T + O2P2.T = 0.
 554   On se place dans le repere de la parabole; soit:
 555          P2 (y*y/(2*p), y, 0)
 556          T (y/p, 1, 0)
 557          D (Dx,Dy,Dz),
 558          V (Vx,Vy,Vz) = (O2O1.D)D-O2O1;
 559
 560   Alors, on obtient l'equation en y suivante:
 561      ((1-Dx*Dx)/(2*p*p))            *  y*y*y  +        A1
 562      (-3*Dx*Dy/(2*p))               *  y*y    +        A2
 563      (1-Dy*Dy + Vx/p)               *  y      +        A3
 564         Vy                          = 0.               A4
 565
 566      On utilise l'algorithme math_DirectPolynomialRoots pour resoudre
 567     cette equation.
 568 -----------------------------------------------------------------------------*/
 569   myIsPar = Standard_False;
 570   myDone = Standard_False;
 571   myNbExt = 0;
 572
 573 // Calcul de T1 dans le repere de la parabole ...
 574   gp_Dir D = C1.Direction();
 575   gp_Dir D1 = D;
 576   gp_Dir x2, y2, z2;
 577   x2 = C2.XAxis().Direction();
 578   y2 = C2.YAxis().Direction();
 579   z2 = C2.Axis().Direction();
 580   Standard_Real Dx = D.Dot(x2);
 581   Standard_Real Dy = D.Dot(y2);
 582   Standard_Real Dz = D.Dot(z2);
 583   D.SetCoord(Dx,Dy,Dz);
 584
 585 // Calcul de V ...
 586   gp_Pnt O1 = C1.Location();
 587   gp_Pnt O2 = C2.Location();
 588   gp_Vec O2O1 (O2,O1);
 589   O2O1.SetCoord(O2O1.Dot(x2), O2O1.Dot(y2), O2O1.Dot(z2));
 590   gp_XYZ Vxyz = (D.XYZ()*(O2O1.Dot(D)))-O2O1.XYZ();
 591
 592 // Calcul des coefficients de l equation en y
 593   Standard_Real P = C2.Parameter();
 594   Standard_Real A1 = (1-Dx*Dx)/(2.0*P*P);
 595   Standard_Real A2 = (-3.0*Dx*Dy/(2.0*P));
 596   Standard_Real A3 = (1 - Dy*Dy + Vxyz.X()/P);
 597   Standard_Real A4 = Vxyz.Y();
 598
 599   math_DirectPolynomialRoots Sol(A1,A2,A3,A4);
 600   if (!Sol.IsDone()) { return; }
 601
 602 // Stockage des solutions ...
 603   gp_Pnt P1,P2;
 604   Standard_Real U1,U2;
 605   Standard_Integer NbSol = Sol.NbSolutions();
 606   for (Standard_Integer NoSol = 1; NoSol <= NbSol; NoSol++) {
 607     U2 = Sol.Value(NoSol);
 608     P2 = ElCLib::Value(U2,C2);
 609     U1 = (gp_Vec(O1,P2)).Dot(D1);
 610     P1 = ElCLib::Value(U1,C1);
 611     mySqDist[myNbExt] = P1.SquareDistance(P2);
 612     myPoint[myNbExt][0] = Extrema_POnCurv(U1,P1);
 613     myPoint[myNbExt][1] = Extrema_POnCurv(U2,P2);
 614     myNbExt++;
 615   }
 616   myDone = Standard_True;
 617 }
 618 //=============================================================================
 619 Extrema_ExtElC::Extrema_ExtElC (const gp_Circ&, const gp_Elips&)
 620 {
 621   Standard_NotImplemented::Raise();
 622 }
 623 //=============================================================================
 624
 625 Extrema_ExtElC::Extrema_ExtElC (const gp_Circ&, const gp_Hypr&)
 626 {
 627   Standard_NotImplemented::Raise();
 628 }
 629 //=============================================================================
 630
 631 Extrema_ExtElC::Extrema_ExtElC (const gp_Circ&, const gp_Parab&)
 632 {
 633   Standard_NotImplemented::Raise();
 634 }
 635 //=============================================================================
 636
 637 Extrema_ExtElC::Extrema_ExtElC (const gp_Elips&, const gp_Elips&)
 638 {
 639   Standard_NotImplemented::Raise();
 640 }
 641 //=============================================================================
 642
 643 Extrema_ExtElC::Extrema_ExtElC (const gp_Elips&, const gp_Hypr&)
 644 {
 645   Standard_NotImplemented::Raise();
 646 }
 647 //=============================================================================
 648
 649 Extrema_ExtElC::Extrema_ExtElC (const gp_Elips&, const gp_Parab&)
 650 {
 651   Standard_NotImplemented::Raise();
 652 }
 653 //=============================================================================
 654
 655 Extrema_ExtElC::Extrema_ExtElC (const gp_Hypr&, const gp_Hypr&)
 656 {
 657   Standard_NotImplemented::Raise();
 658 }
 659 //=============================================================================
 660
 661 Extrema_ExtElC::Extrema_ExtElC (const gp_Hypr&, const gp_Parab&)
 662 {
 663   Standard_NotImplemented::Raise();
 664 }
 665 //=============================================================================
 666
 667 Extrema_ExtElC::Extrema_ExtElC (const gp_Parab&, const gp_Parab&)
 668 {
 669   Standard_NotImplemented::Raise();
 670 }
 671 //=============================================================================
 672
 673 Standard_Boolean Extrema_ExtElC::IsDone () const { return myDone; }
 674 //=============================================================================
 675
 676 Standard_Boolean Extrema_ExtElC::IsParallel () const
 677 {
 678   if (!IsDone()) { StdFail_NotDone::Raise(); }
 679   return myIsPar;
 680 }
 681 //=============================================================================
 682
 683 Standard_Integer Extrema_ExtElC::NbExt () const
 684 {
 685   if (IsParallel()) { StdFail_InfiniteSolutions::Raise(); }
 686   return myNbExt;
 687 }
 688 //=============================================================================
 689
 690 Standard_Real Extrema_ExtElC::SquareDistance (const Standard_Integer N) const
 691 {
 692   if (!myDone) { StdFail_NotDone::Raise(); }
 693   if (myIsPar) {
 694     if (N < 1 || N > 2) { Standard_OutOfRange::Raise(); }
 695   }
 696   else {
 697     if (N < 1 || N > NbExt()) { Standard_OutOfRange::Raise(); }
 698   }
 699   return mySqDist[N-1];
 700 }
 701 //=============================================================================
 702
 703 void Extrema_ExtElC::Points (const Standard_Integer N,
 704                              Extrema_POnCurv& P1, Extrema_POnCurv& P2) const
 705 {
 706   if (N < 1 || N > NbExt()) { Standard_OutOfRange::Raise(); }
 707   P1 = myPoint[N-1][0];
 708   P2 = myPoint[N-1][1];
 709 }
 710 //=============================================================================
 711 //=============================================================================
 712 //
 713 //modified by NIZNHY-PKV Fri Nov 21 10:48:46 2008f
 714 //=======================================================================
 715 //function : Extrema_ExtElC
 716 //purpose  :
 717 //=======================================================================
 718 Extrema_ExtElC::Extrema_ExtElC (const gp_Circ& C1, const gp_Circ& C2)
 719 {
 720   Standard_Boolean bIsSamePlane, bIsSameAxe;
 721   Standard_Real aTolD, aTolD2, aTolA, aD2, aDC2;
 722   gp_Pnt aPc1, aPc2;
 723   gp_Dir aDc1, aDc2;
 724   //
 725   myIsPar = Standard_False;
 726   myDone = Standard_False;
 727   myNbExt = 0;
 728   //
 729   aTolA=Precision::Angular();
 730   aTolD=Precision::Confusion();
 731   aTolD2=aTolD*aTolD;
 732   //
 733   aPc1=C1.Location();
 734   aDc1=C1.Axis().Direction();
 735
 736   aPc2=C2.Location();
 737   aDc2=C2.Axis().Direction();
 738   gp_Pln aPlc1(aPc1, aDc1);
 739   //
 740   aD2=aPlc1.SquareDistance(aPc2);
 741   bIsSamePlane=aDc1.IsParallel(aDc2, aTolA) && aD2<aTolD2;
 742   if (!bIsSamePlane) {
 743     return;
 744   }
 745   //
 746   aDC2=aPc1.SquareDistance(aPc2);
 747   bIsSameAxe=aDC2<aTolD2;
 748   //
 749   if(bIsSameAxe) {
 750     myIsPar = Standard_True;
 751     Standard_Real dR = C1.Radius() - C2.Radius();
 752     Standard_Real dC = C1.Location().Distance(C2.Location());
 753     mySqDist[0] = dR*dR + dC*dC;
 754     dR = C1.Radius() + C2.Radius();
 755     mySqDist[1] = dR*dR + dC*dC;
 756     myDone = Standard_True;
 757   }
 758   else {
 759     Standard_Boolean bIn, bOut;
 760     Standard_Integer j1, j2;
 761     Standard_Real aR1, aR2, aD12, aT11, aT12, aT21, aT22;
 762     gp_Circ aC1, aC2;
 763     gp_Pnt aP11, aP12, aP21, aP22;
 764     //
 765     myDone = Standard_True;
 766     //
 767     aR1=C1.Radius();
 768     aR2=C2.Radius();
 769     //
 770     j1=0;
 771     j2=1;
 772     aC1=C1;
 773     aC2=C2;
 774     if (aR2>aR1) {
 775       j1=1;
 776       j2=0;
 777       aC1=C2;
 778       aC2=C1;
 779     }
 780     //
 781     aR1=aC1.Radius(); // max radius
 782     aR2=aC2.Radius(); // min radius
 783     //
 784     aPc1=aC1.Location();
 785     aPc2=aC2.Location();
 786     //
 787     aD12=aPc1.Distance(aPc2);
 788     gp_Vec aVec12(aPc1, aPc2);
 789     gp_Dir aDir12(aVec12);
 790     //
 791     // 1. Four common solutions
 792     myNbExt=4;
 793     //
 794     aP11.SetXYZ(aPc1.XYZ()-aR1*aDir12.XYZ());
 795     aP12.SetXYZ(aPc1.XYZ()+aR1*aDir12.XYZ());
 796     aP21.SetXYZ(aPc2.XYZ()-aR2*aDir12.XYZ());
 797     aP22.SetXYZ(aPc2.XYZ()+aR2*aDir12.XYZ());
 798     //
 799     aT11=ElCLib::Parameter(aC1, aP11);
 800     aT12=ElCLib::Parameter(aC1, aP12);
 801     aT21=ElCLib::Parameter(aC2, aP21);
 802     aT22=ElCLib::Parameter(aC2, aP22);
 803     //
 804     // P11, P21
 805     myPoint[0][j1].SetValues(aT11, aP11);
 806     myPoint[0][j2].SetValues(aT21, aP21);
 807     mySqDist[0]=aP11.SquareDistance(aP21);
 808     // P11, P22
 809     myPoint[1][j1].SetValues(aT11, aP11);
 810     myPoint[1][j2].SetValues(aT22, aP22);
 811     mySqDist[1]=aP11.SquareDistance(aP22);
 812     //
 813     // P12, P21
 814     myPoint[2][j1].SetValues(aT12, aP12);
 815     myPoint[2][j2].SetValues(aT21, aP21);
 816     mySqDist[2]=aP12.SquareDistance(aP21);
 817     //
 818     // P12, P22
 819     myPoint[3][j1].SetValues(aT12, aP12);
 820     myPoint[3][j2].SetValues(aT22, aP22);
 821     mySqDist[3]=aP12.SquareDistance(aP22);
 822     //
 823     // 2. Check for intersections
 824     bOut=aD12>(aR1+aR2+aTolD);
 825     bIn =aD12<(aR1-aR2-aTolD);
 826     if (!bOut && !bIn) {
 827       Standard_Boolean bNbExt6;
 828       Standard_Real aAlpha, aBeta, aT[2], aVal, aDist2;
 829       gp_Pnt aPt, aPL1, aPL2;
 830       gp_Dir aDLt;
 831       //
 832       aAlpha=0.5*(aR1*aR1-aR2*aR2+aD12*aD12)/aD12;
 833       aVal=aR1*aR1-aAlpha*aAlpha;
 834       if (aVal<0.) {// see pkv/900/L4 for details
 835         aVal=-aVal;
 836       }
 837       aBeta=Sqrt(aVal);
 838       //aBeta=Sqrt(aR1*aR1-aAlpha*aAlpha);
 839       //--
 840       aPt.SetXYZ(aPc1.XYZ()+aAlpha*aDir12.XYZ());
 841       //
 842       aDLt=aDc1^aDir12;
 843       aPL1.SetXYZ(aPt.XYZ()+aBeta*aDLt.XYZ());
 844       aPL2.SetXYZ(aPt.XYZ()-aBeta*aDLt.XYZ());
 845       //
 846       aDist2=aPL1.SquareDistance(aPL2);
 847       bNbExt6=aDist2>aTolD2;
 848       //
 849       myNbExt=5;// just in case. see pkv/900/L4 for details
 850       aT[j1]=ElCLib::Parameter(aC1, aPL1);
 851       aT[j2]=ElCLib::Parameter(aC2, aPL1);
 852       myPoint[4][j1].SetValues(aT[j1], aPL1);
 853       myPoint[4][j2].SetValues(aT[j2], aPL1);
 854       mySqDist[4]=0.;
 855       //
 856       if (bNbExt6) {
 857         myNbExt=6;
 858         aT[j1]=ElCLib::Parameter(aC1, aPL2);
 859         aT[j2]=ElCLib::Parameter(aC2, aPL2);
 860         myPoint[5][j1].SetValues(aT[j1], aPL2);
 861         myPoint[5][j2].SetValues(aT[j2], aPL2);
 862         mySqDist[5]=0.;
 863       }
 864       //
 865     }// if (!bOut || !bIn) {
 866   }// else
 867 }
 868 //modified by NIZNHY-PKV Fri Nov 21 10:48:56 2008t