2fcbc1e3de3a216feef979f1e3b9f6e9cd3f033e
[occt.git] / 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