0024915: Wrong intersection curves between two cylinders
[occt.git] / src / IntAna / IntAna_QuadQuadGeo.cxx
1 // Created on: 1992-08-06
2 // Created by: Laurent BUCHARD
3 // Copyright (c) 1992-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
8 // This library is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU Lesser General Public License version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17 //----------------------------------------------------------------------
18 //-- Purposse: Geometric Intersection between two Natural Quadric 
19 //--          If the intersection is not a conic, 
20 //--          analytical methods must be called.
21 //----------------------------------------------------------------------
22 #ifndef DEB
23 #define No_Standard_RangeError
24 #define No_Standard_OutOfRange
25 #endif
26
27 #include <IntAna_QuadQuadGeo.ixx>
28
29 #include <IntAna_IntConicQuad.hxx>
30 #include <StdFail_NotDone.hxx>
31 #include <Standard_DomainError.hxx>
32 #include <Standard_OutOfRange.hxx>
33 #include <math_DirectPolynomialRoots.hxx>
34
35 #include <gp.hxx>
36 #include <gp_Pln.hxx>
37 #include <gp_Vec.hxx>
38 #include <ElSLib.hxx>
39 #include <ElCLib.hxx>
40
41 #include <gp_Dir.hxx>
42 #include <gp_XYZ.hxx>
43 #include <gp_Pnt2d.hxx>
44 #include <gp_Vec2d.hxx>
45 #include <gp_Dir2d.hxx>
46
47
48 static
49   gp_Ax2 DirToAx2(const gp_Pnt& P,const gp_Dir& D);
50 static
51   void RefineDir(gp_Dir& aDir);
52
53 //=======================================================================
54 //class :  AxeOperator
55 //purpose  : O p e r a t i o n s   D i v e r s e s  s u r   d e s   A x 1 
56 //=======================================================================
57 class AxeOperator {
58  public:
59   AxeOperator(const gp_Ax1& A1,const gp_Ax1& A2);
60
61   void Distance(Standard_Real& dist,
62                 Standard_Real& Param1,
63                 Standard_Real& Param2);
64   
65   gp_Pnt PtIntersect()              { 
66     return ptintersect;
67   }
68   Standard_Boolean Coplanar(void)   { 
69     return thecoplanar;
70   }
71   Standard_Boolean Same(void)       {
72     return theparallel && (thedistance<myEPSILON_DISTANCE); 
73   }
74   Standard_Real Distance(void)      { 
75     return thedistance ;
76   }
77   Standard_Boolean Intersect(void)  { 
78     return (thecoplanar && (!theparallel));
79   }
80   Standard_Boolean Parallel(void)   { 
81     return theparallel; 
82   }
83   Standard_Boolean Normal(void)     { 
84     return thenormal;
85   }
86
87  protected:
88   Standard_Real Det33(const Standard_Real a11,
89                       const Standard_Real a12,
90                       const Standard_Real a13,
91                       const Standard_Real a21,
92                       const Standard_Real a22,
93                       const Standard_Real a23,
94                       const Standard_Real a31,
95                       const Standard_Real a32,
96                       const Standard_Real a33) {
97     Standard_Real theReturn =  
98       a11*(a22*a33-a32*a23) - a21*(a12*a33-a32*a13) + a31*(a12*a23-a22*a13) ;   
99     return theReturn ;
100   }
101
102  private:
103   gp_Pnt ptintersect;
104   gp_Ax1 Axe1;
105   gp_Ax1 Axe2;
106   Standard_Real thedistance;
107   Standard_Boolean theparallel;
108   Standard_Boolean thecoplanar;
109   Standard_Boolean thenormal;
110   //
111   Standard_Real myEPSILON_DISTANCE;
112   Standard_Real myEPSILON_AXES_PARA;
113 };
114
115 //=======================================================================
116 //function : AxeOperator::AxeOperator
117 //purpose  : 
118 //=======================================================================
119 AxeOperator::AxeOperator(const gp_Ax1& A1,const gp_Ax1& A2) 
120 {
121   myEPSILON_DISTANCE=0.00000000000001;
122   myEPSILON_AXES_PARA=0.000000000001;
123   Axe1=A1; 
124   Axe2=A2;
125   //---------------------------------------------------------------------
126   gp_Dir V1=Axe1.Direction();
127   gp_Dir V2=Axe2.Direction();
128   gp_Pnt P1=Axe1.Location();
129   gp_Pnt P2=Axe2.Location();
130   //
131   RefineDir(V1);
132   RefineDir(V2);
133   thecoplanar= Standard_False;
134   thenormal  = Standard_False;
135
136   //--- check if the two axis are parallel
137   theparallel=V1.IsParallel(V2, myEPSILON_AXES_PARA);  
138   //--- Distance between the two axis
139   gp_XYZ perp(A1.Direction().XYZ().Crossed(A2.Direction().XYZ()));
140   if (theparallel) { 
141     gp_Lin L1(A1);
142     thedistance = L1.Distance(A2.Location());
143   }
144   else {   
145     thedistance = Abs(gp_Vec(perp.Normalized()).Dot(gp_Vec(Axe1.Location(),
146                                                            Axe2.Location())));
147   }
148   //--- check if Axis are Coplanar
149   Standard_Real D33;
150   if(thedistance<myEPSILON_DISTANCE) {
151     D33=Det33(V1.X(),V1.Y(),V1.Z()
152               ,V2.X(),V2.Y(),V2.Z()
153               ,P1.X()-P2.X(),P1.Y()-P2.Y(),P1.Z()-P2.Z());
154     if(Abs(D33)<=myEPSILON_DISTANCE) { 
155       thecoplanar=Standard_True;
156     }
157   }
158   else {
159     thecoplanar=Standard_True;
160     thenormal=(V1.Dot(V2)==0.0)? Standard_True : Standard_False;
161   }
162   //--- check if the two axis are concurrent
163   if(thecoplanar && (!theparallel)) {
164     Standard_Real smx=P2.X() - P1.X();
165     Standard_Real smy=P2.Y() - P1.Y();
166     Standard_Real smz=P2.Z() - P1.Z();
167     Standard_Real Det1,Det2,Det3,A;
168     Det1=V1.Y() * V2.X() - V1.X() * V2.Y();
169     Det2=V1.Z() * V2.Y() - V1.Y() * V2.Z();
170     Det3=V1.Z() * V2.X() - V1.X() * V2.Z();
171     
172     if((Det1!=0.0) && ((Abs(Det1) >= Abs(Det2))&&(Abs(Det1) >= Abs(Det3)))) {
173       A=(smy * V2.X() - smx * V2.Y())/Det1;
174     }
175     else if((Det2!=0.0) 
176             && ((Abs(Det2) >= Abs(Det1))
177                 &&(Abs(Det2) >= Abs(Det3)))) {
178       A=(smz * V2.Y() - smy * V2.Z())/Det2;
179     }
180     else {
181       A=(smz * V2.X() - smx * V2.Z())/Det3;
182     }
183     ptintersect.SetCoord( P1.X() + A * V1.X()
184                          ,P1.Y() + A * V1.Y()
185                          ,P1.Z() + A * V1.Z());
186   }
187   else { 
188     ptintersect.SetCoord(0,0,0);  //-- Pour eviter des FPE
189   }
190 }
191 //=======================================================================
192 //function : Distance
193 //purpose  : 
194 //=======================================================================
195 void AxeOperator::Distance(Standard_Real& dist,
196                            Standard_Real& Param1,
197                            Standard_Real& Param2)
198  {
199   gp_Vec O1O2(Axe1.Location(),Axe2.Location());   
200   gp_Dir U1 = Axe1.Direction();   //-- juste pour voir. 
201   gp_Dir U2 = Axe2.Direction();
202   
203   gp_Dir N  = U1.Crossed(U2);
204   Standard_Real D = Det33(U1.X(),U2.X(),N.X(),
205                           U1.Y(),U2.Y(),N.Y(),
206                           U1.Z(),U2.Z(),N.Z());
207   if(D) { 
208     dist = Det33(U1.X(),U2.X(),O1O2.X(),
209                  U1.Y(),U2.Y(),O1O2.Y(),
210                  U1.Z(),U2.Z(),O1O2.Z()) / D;
211     Param1 = Det33(O1O2.X(),U2.X(),N.X(),
212                    O1O2.Y(),U2.Y(),N.Y(),
213                    O1O2.Z(),U2.Z(),N.Z()) / (-D);
214     //------------------------------------------------------------
215     //-- On resout P1 * Dir1 + P2 * Dir2 + d * N = O1O2
216     //-- soit : Segment perpendiculaire : O1+P1 D1
217     //--                                  O2-P2 D2
218     Param2 = Det33(U1.X(),O1O2.X(),N.X(),
219                    U1.Y(),O1O2.Y(),N.Y(),
220                    U1.Z(),O1O2.Z(),N.Z()) / (D);
221   }
222 }
223 //=======================================================================
224 //function : DirToAx2
225 //purpose  : returns a gp_Ax2 where D is the main direction
226 //=======================================================================
227 gp_Ax2 DirToAx2(const gp_Pnt& P,const gp_Dir& D) 
228 {
229   Standard_Real x=D.X(); Standard_Real ax=Abs(x);
230   Standard_Real y=D.Y(); Standard_Real ay=Abs(y);
231   Standard_Real z=D.Z(); Standard_Real az=Abs(z);
232   if( (ax==0.0) || ((ax<ay) && (ax<az)) ) {
233     return(gp_Ax2(P,D,gp_Dir(gp_Vec(0.0,-z,y))));
234   }
235   else if( (ay==0.0) || ((ay<ax) && (ay<az)) ) {
236     return(gp_Ax2(P,D,gp_Dir(gp_Vec(-z,0.0,x))));
237   }
238   else {
239     return(gp_Ax2(P,D,gp_Dir(gp_Vec(-y,x,0.0))));
240   }
241 }
242 //=======================================================================
243 //function : IntAna_QuadQuadGeo
244 //purpose  : Empty constructor
245 //=======================================================================
246 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(void)
247     : done(Standard_False),
248       nbint(0),
249       typeres(IntAna_Empty),
250       pt1(0,0,0),
251       pt2(0,0,0),
252       pt3(0,0,0),
253       pt4(0,0,0),
254       param1(0),
255       param2(0),
256       param3(0),
257       param4(0),
258       param1bis(0),
259       param2bis(0),
260       myCommonGen(Standard_False),
261       myPChar(0,0,0)
262 {
263   InitTolerances();
264 }
265 //=======================================================================
266 //function : InitTolerances
267 //purpose  : 
268 //=======================================================================
269 void IntAna_QuadQuadGeo::InitTolerances()
270 {
271   myEPSILON_DISTANCE               = 0.00000000000001;
272   myEPSILON_ANGLE_CONE             = 0.000000000001;
273   myEPSILON_MINI_CIRCLE_RADIUS     = 0.000000001;
274   myEPSILON_CYLINDER_DELTA_RADIUS  = 0.0000000000001;
275   myEPSILON_CYLINDER_DELTA_DISTANCE= 0.0000001;
276   myEPSILON_AXES_PARA              = 0.000000000001;
277 }
278 //=======================================================================
279 //function : IntAna_QuadQuadGeo
280 //purpose  : Pln  Pln 
281 //=======================================================================
282 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P1, 
283                                        const gp_Pln& P2,
284                                        const Standard_Real TolAng,
285                                        const Standard_Real Tol)
286 : done(Standard_False),
287   nbint(0),
288   typeres(IntAna_Empty),
289   pt1(0,0,0),
290   pt2(0,0,0),
291   pt3(0,0,0),
292   pt4(0,0,0),
293   param1(0),
294   param2(0),
295   param3(0),
296   param4(0),
297   param1bis(0),
298   param2bis(0),
299   myCommonGen(Standard_False),
300   myPChar(0,0,0)
301 {
302   InitTolerances();
303   Perform(P1,P2,TolAng,Tol);
304 }
305 //=======================================================================
306 //function : Perform
307 //purpose  : 
308 //=======================================================================
309 void IntAna_QuadQuadGeo::Perform (const gp_Pln& P1, 
310                                   const gp_Pln& P2,
311                                   const Standard_Real TolAng,
312                                   const Standard_Real Tol)
313 {
314   Standard_Real A1, B1, C1, D1, A2, B2, C2, D2, dist1, dist2, aMVD;
315   //
316   done=Standard_False;
317   param2bis=0.;
318   //
319   P1.Coefficients(A1,B1,C1,D1);
320   P2.Coefficients(A2,B2,C2,D2);
321   //
322   gp_Vec aVN1(A1,B1,C1);
323   gp_Vec aVN2(A2,B2,C2);
324   gp_Vec vd(aVN1.Crossed(aVN2));
325   //
326   const gp_Pnt& aLocP1=P1.Location();
327   const gp_Pnt& aLocP2=P2.Location();
328   //
329   dist1=A2*aLocP1.X() + B2*aLocP1.Y() + C2*aLocP1.Z() + D2;
330   dist2=A1*aLocP2.X() + B1*aLocP2.Y() + C1*aLocP2.Z() + D1;
331   //
332   aMVD=vd.Magnitude();
333   if(aMVD <=TolAng) {
334     // normalles are collinear - planes are same or parallel
335     typeres = (Abs(dist1) <= Tol && Abs(dist2) <= Tol) ? IntAna_Same 
336       : IntAna_Empty;
337   }
338   else {
339     Standard_Real denom, denom2, ddenom, par1, par2;
340     Standard_Real X1, Y1, Z1, X2, Y2, Z2, aEps;
341     //
342     aEps=1.e-16;
343     denom=A1*A2 + B1*B2 + C1*C2;
344     denom2 = denom*denom;
345     ddenom = 1. - denom2;
346     
347     denom = ( Abs(ddenom) <= aEps ) ? aEps : ddenom;
348       
349     par1 = dist1/denom;
350     par2 = -dist2/denom;
351       
352     gp_Vec inter1(aVN1.Crossed(vd));
353     gp_Vec inter2(aVN2.Crossed(vd));
354       
355     X1=aLocP1.X() + par1*inter1.X();
356     Y1=aLocP1.Y() + par1*inter1.Y();
357     Z1=aLocP1.Z() + par1*inter1.Z();
358     X2=aLocP2.X() + par2*inter2.X();
359     Y2=aLocP2.Y() + par2*inter2.Y();
360     Z2=aLocP2.Z() + par2*inter2.Z();
361       
362     pt1=gp_Pnt((X1+X2)*0.5, (Y1+Y2)*0.5, (Z1+Z2)*0.5);
363     dir1 = gp_Dir(vd);
364     typeres = IntAna_Line;
365     nbint = 1;
366     //
367     //-------------------------------------------------------
368     // When the value of the angle between the planes is small
369     // the origin of intersection line is computed with error
370     // [ ~0.0001 ] that can not br considered as small one
371     // e.g.
372     // for {A~=2.e-6, dist1=4.2e-5, dist2==1.e-4} =>
373     // {denom=3.4e-12, par1=12550297.6, par2=32605552.9, etc}
374     // So, 
375     // the origin should be refined if it is possible
376     //
377     Standard_Real aTreshAng, aTreshDist;
378     //
379     aTreshAng=2.e-6; // 1.e-4 deg
380     aTreshDist=1.e-12;
381     //
382     if (aMVD < aTreshAng) {
383       Standard_Real aDist1, aDist2;
384       //
385       aDist1=A1*pt1.X() + B1*pt1.Y() + C1*pt1.Z() + D1;
386       aDist2=A2*pt1.X() + B2*pt1.Y() + C2*pt1.Z() + D2;
387       //
388       if (fabs(aDist1)>aTreshDist || fabs(aDist2)>aTreshDist) {
389         Standard_Boolean bIsDone, bIsParallel;
390         IntAna_IntConicQuad aICQ;
391         //
392         // 1.
393         gp_Dir aDN1(aVN1);
394         gp_Lin aL1(pt1, aDN1);
395         //
396         aICQ.Perform(aL1, P1, TolAng, Tol);
397         bIsDone=aICQ.IsDone();
398         if (!bIsDone) {
399           return;
400         }
401         //
402         const gp_Pnt& aPnt1=aICQ.Point(1);
403         //----------------------------------
404         // 2.
405         gp_Dir aDL2(dir1.Crossed(aDN1));
406         gp_Lin aL2(aPnt1, aDL2);
407         //
408         aICQ.Perform(aL2, P2, TolAng, Tol);
409         bIsDone=aICQ.IsDone();
410         if (!bIsDone) {
411           return;
412         }
413         //
414         bIsParallel=aICQ.IsParallel();
415         if (bIsParallel) {
416           return;
417         }
418         //
419         const gp_Pnt& aPnt2=aICQ.Point(1);
420         //
421         pt1=aPnt2;
422       }
423     }
424   }
425   done=Standard_True;
426 }
427 //=======================================================================
428 //function : IntAna_QuadQuadGeo
429 //purpose  : Pln Cylinder
430 //=======================================================================
431 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo( const gp_Pln& P
432                                        ,const gp_Cylinder& Cl
433                                        ,const Standard_Real Tolang
434                                        ,const Standard_Real Tol
435                                        ,const Standard_Real H)
436      : done(Standard_False),
437        nbint(0),
438        typeres(IntAna_Empty),
439        pt1(0,0,0),
440        pt2(0,0,0),
441        pt3(0,0,0),
442        pt4(0,0,0),
443        param1(0),
444        param2(0),
445        param3(0),
446        param4(0),
447        param1bis(0),
448        param2bis(0),
449        myCommonGen(Standard_False),
450        myPChar(0,0,0)
451 {
452   InitTolerances();
453   Perform(P,Cl,Tolang,Tol,H);
454 }
455 //=======================================================================
456 //function : Perform
457 //purpose  : 
458 //=======================================================================
459   void IntAna_QuadQuadGeo::Perform( const gp_Pln& P
460                                    ,const gp_Cylinder& Cl
461                                    ,const Standard_Real Tolang
462                                    ,const Standard_Real Tol
463                                    ,const Standard_Real H) 
464 {
465   done = Standard_False;
466   Standard_Real dist,radius;
467   Standard_Real A,B,C,D;
468   Standard_Real X,Y,Z;
469   Standard_Real sint,cost,h;
470   gp_XYZ axex,axey,omega;
471
472   
473   param2bis=0.0;
474   radius = Cl.Radius();
475
476   gp_Lin axec(Cl.Axis());
477   gp_XYZ normp(P.Axis().Direction().XYZ());
478
479   P.Coefficients(A,B,C,D);
480   axec.Location().Coord(X,Y,Z);
481   // la distance axe/plan est evaluee a l origine de l axe.
482   dist = A*X + B*Y + C*Z + D; 
483
484   Standard_Real tolang = Tolang;
485   Standard_Boolean newparams = Standard_False;
486
487   gp_Vec ldv( axec.Direction() );
488   gp_Vec npv( normp );
489   Standard_Real dA = Abs( ldv.Angle( npv ) );
490   if( dA > (M_PI/4.) )
491     {
492       Standard_Real dang = Abs( ldv.Angle( npv ) ) - M_PI/2.;
493       Standard_Real dangle = Abs( dang );
494       if( dangle > Tolang )
495         {
496           Standard_Real sinda = Abs( Sin( dangle ) );
497           Standard_Real dif = Abs( sinda - Tol );
498           if( dif < Tol )
499             {
500               tolang = sinda * 2.;
501               newparams = Standard_True;
502             }
503         }
504     }
505
506   nbint = 0;
507   IntAna_IntConicQuad inter(axec,P,tolang,Tol,H);
508
509   if (inter.IsParallel()) {
510     // Le resultat de l intersection Plan-Cylindre est de type droite.
511     // il y a 1 ou 2 droites
512
513     typeres = IntAna_Line;
514     omega.SetCoord(X-dist*A,Y-dist*B,Z-dist*C);
515
516     if (Abs(Abs(dist)-radius) < Tol)
517       {
518         nbint = 1;
519         pt1.SetXYZ(omega);
520
521         if( newparams )
522           {
523             gp_XYZ omegaXYZ(X,Y,Z);
524             gp_XYZ omegaXYZtrnsl( omegaXYZ + 100.*axec.Direction().XYZ() );
525             Standard_Real Xt,Yt,Zt,distt;
526             omegaXYZtrnsl.Coord(Xt,Yt,Zt);
527             distt = A*Xt + B*Yt + C*Zt + D;
528             gp_XYZ omega1(omegaXYZtrnsl.X()-distt*A, 
529                           omegaXYZtrnsl.Y()-distt*B, 
530                           omegaXYZtrnsl.Z()-distt*C );
531             gp_Pnt ppt1;
532             ppt1.SetXYZ( omega1 );
533             gp_Vec vv1(pt1,ppt1);
534             gp_Dir dd1( vv1 );
535             dir1 = dd1;
536           }
537         else
538           dir1 = axec.Direction();
539     }
540     else if (Abs(dist) < radius)
541       {
542         nbint = 2;
543         h = Sqrt(radius*radius - dist*dist);
544         axey = axec.Direction().XYZ().Crossed(normp); // axey est normalise
545
546         pt1.SetXYZ(omega - h*axey);
547         pt2.SetXYZ(omega + h*axey);
548
549         if( newparams )
550           { 
551             gp_XYZ omegaXYZ(X,Y,Z);
552             gp_XYZ omegaXYZtrnsl( omegaXYZ + 100.*axec.Direction().XYZ() );
553             Standard_Real Xt,Yt,Zt,distt,ht;
554             omegaXYZtrnsl.Coord(Xt,Yt,Zt);
555             distt = A*Xt + B*Yt + C*Zt + D;
556             //             ht = Sqrt(radius*radius - distt*distt);
557             Standard_Real anSqrtArg = radius*radius - distt*distt;
558             ht = (anSqrtArg > 0.) ? Sqrt(anSqrtArg) : 0.;
559
560             gp_XYZ omega1( omegaXYZtrnsl.X()-distt*A, 
561                           omegaXYZtrnsl.Y()-distt*B, 
562                           omegaXYZtrnsl.Z()-distt*C );
563             gp_Pnt ppt1,ppt2;
564             ppt1.SetXYZ( omega1 - ht*axey);
565             ppt2.SetXYZ( omega1 + ht*axey);
566             gp_Vec vv1(pt1,ppt1);
567             gp_Vec vv2(pt2,ppt2);
568             gp_Dir dd1( vv1 );
569             gp_Dir dd2( vv2 );
570             dir1 = dd1;
571             dir2 = dd2;
572           }
573         else
574           {
575             dir1 = axec.Direction();
576             dir2 = axec.Direction();
577           }
578     }
579     //  else nbint = 0
580
581     // debug JAG : le nbint = 0 doit etre remplace par typeres = IntAna_Empty
582     // et ne pas etre seulement supprime...
583
584     else {
585       typeres = IntAna_Empty;
586     }
587   }
588   else {     // Il y a un point d intersection. C est le centre du cercle
589              // ou de l ellipse solution.
590
591     nbint = 1;
592     axey = normp.Crossed(axec.Direction().XYZ());
593     sint = axey.Modulus();
594
595     pt1 = inter.Point(1);
596     
597     if (sint < Tol/radius) {
598
599       // on construit un cercle avec comme axes X et Y ceux du cylindre
600       typeres = IntAna_Circle;
601
602       dir1 = axec.Direction(); // axe Z
603       dir2 = Cl.Position().XDirection();
604       param1 = radius;
605     }
606     else {
607
608       // on construit un ellipse
609       typeres = IntAna_Ellipse;
610       cost = Abs(axec.Direction().XYZ().Dot(normp));
611       axex = axey.Crossed(normp);
612
613       dir1.SetXYZ(normp);         //Modif ds ce bloc 
614       dir2.SetXYZ(axex);
615
616       param1    = radius/cost;
617       param1bis = radius;
618     }
619   }
620
621   done = Standard_True;
622 }
623 //=======================================================================
624 //function : IntAna_QuadQuadGeo
625 //purpose  : Pln Cone
626 //=======================================================================
627   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P,
628                                          const gp_Cone& Co,
629                                          const Standard_Real Tolang,
630                                          const Standard_Real Tol) 
631 : done(Standard_False),
632   nbint(0),
633   typeres(IntAna_Empty),
634   pt1(0,0,0),
635   pt2(0,0,0),
636   pt3(0,0,0),
637   pt4(0,0,0),
638   param1(0),
639   param2(0),
640   param3(0),
641   param4(0),
642   param1bis(0),
643   param2bis(0),
644   myCommonGen(Standard_False),
645   myPChar(0,0,0)
646 {
647   InitTolerances();
648   Perform(P,Co,Tolang,Tol);
649 }
650 //=======================================================================
651 //function : Perform
652 //purpose  : 
653 //=======================================================================
654   void IntAna_QuadQuadGeo::Perform(const gp_Pln& P,
655                                    const gp_Cone& Co,
656                                    const Standard_Real Tolang,
657                                    const Standard_Real Tol) 
658 {
659
660   done = Standard_False;
661   nbint = 0;
662
663   Standard_Real A,B,C,D;
664   Standard_Real X,Y,Z;
665   Standard_Real dist,sint,cost,sina,cosa,angl,costa;
666   Standard_Real dh;
667   gp_XYZ axex,axey;
668
669   gp_Lin axec(Co.Axis());
670   P.Coefficients(A,B,C,D);
671   gp_Pnt apex(Co.Apex());
672
673   apex.Coord(X,Y,Z);
674   dist = A*X + B*Y + C*Z + D; // distance signee sommet du cone/ Plan
675
676   gp_XYZ normp = P.Axis().Direction().XYZ();
677   if(P.Direct()==Standard_False) {  //-- lbr le 14 jan 97
678     normp.Reverse();
679   }
680
681   axey = normp.Crossed(Co.Axis().Direction().XYZ());
682   axex = axey.Crossed(normp);
683
684
685   angl = Co.SemiAngle();
686
687   cosa = Cos(angl);
688   sina = Abs(Sin(angl));
689
690
691   // Angle entre la normale au plan et l axe du cone, ramene entre 0. et PI/2.
692
693   sint = axey.Modulus();
694   cost = Abs(Co.Axis().Direction().XYZ().Dot(normp));
695
696   // Le calcul de costa permet de determiner si le plan contient
697   // un generatrice du cone : on calcul Sin((PI/2. - t) - angl)
698
699   costa = cost*cosa - sint*sina;  // sin((PI/2 -t)-angl)=cos(t+angl)
700                                   // avec  t ramene entre 0 et pi/2.
701
702   if (Abs(dist) < Tol) {
703     // on considere que le plan contient le sommet du cone.
704     // les solutions possibles sont donc : 1 point, 1 droite, 2 droites
705     // selon l inclinaison du plan.
706
707     if (Abs(costa) < Tolang) {          // plan parallele a la generatrice
708       typeres = IntAna_Line;
709       nbint = 1;
710       gp_XYZ ptonaxe(apex.XYZ() + 10.*(Co.Axis().Direction().XYZ()));
711       // point sur l axe du cone cote z positif
712
713       dist = A*ptonaxe.X() + B*ptonaxe.Y() + C*ptonaxe.Z() + D;
714       ptonaxe = ptonaxe - dist*normp;
715       pt1 = apex;
716       dir1.SetXYZ(ptonaxe - pt1.XYZ());
717     }
718     else if (cost < sina) {   // plan "interieur" au cone
719       typeres = IntAna_Line;
720       nbint = 2;
721       pt1 = apex;
722       pt2 = apex;
723       dh = Sqrt(sina*sina-cost*cost)/cosa;
724       dir1.SetXYZ(axex + dh*axey);
725       dir2.SetXYZ(axex - dh*axey);
726     }
727     else {                              // plan "exterieur" au cone
728       typeres = IntAna_Point;
729       nbint = 1;
730       pt1 = apex;
731     }
732   }
733   else {
734     // Solutions possibles : cercle, ellipse, parabole, hyperbole selon
735     // l inclinaison du plan.
736     Standard_Real deltacenter, distance;
737
738     if (cost < Tolang) {
739       // Le plan contient la direction de l axe du cone. La solution est
740       // l hyperbole
741       typeres = IntAna_Hyperbola;
742       nbint = 2;
743       pt1.SetXYZ(apex.XYZ()-dist*normp);
744       pt2 = pt1;
745       dir1=normp;
746       dir2.SetXYZ(axex);
747       param1     = param2 = Abs(dist/Tan(angl));
748       param1bis  = param2bis = Abs(dist);
749     }
750     else {
751
752       IntAna_IntConicQuad inter(axec,P,Tolang); // on a necessairement 1 point.
753       
754       gp_Pnt center(inter.Point(1));
755
756       // En fonction de la position de l intersection par rapport au sommet
757       // du cone, on change l axe x en -x et y en -y. Le parametre du sommet
758       // sur axec est negatif (voir definition du cone)
759
760       distance = apex.Distance(center);
761
762       if (inter.ParamOnConic(1) + Co.RefRadius()/Tan(angl) < 0.) {
763         axex.Reverse();
764         axey.Reverse();
765       }
766
767       if (Abs(costa) < Tolang) {  // plan parallele a une generatrice
768         typeres = IntAna_Parabola;
769         nbint = 1;
770         deltacenter = distance/2./cosa;
771         axex.Normalize();
772         pt1.SetXYZ(center.XYZ()-deltacenter*axex);
773         dir1 = normp;
774         dir2.SetXYZ(axex);
775         param1 = deltacenter*sina*sina;
776       }
777       else if (sint  < Tolang) {            // plan perpendiculaire a l axe
778         typeres = IntAna_Circle;
779         nbint = 1;
780         pt1 = center;
781         dir1 = Co.Position().Direction();
782         dir2 = Co.Position().XDirection();
783         param1 = apex.Distance(center)*Abs(Tan(angl));
784       }
785       else if (cost < sina ) {
786         typeres = IntAna_Hyperbola;
787         nbint = 2;
788         axex.Normalize();
789
790         deltacenter = sint*sina*sina*distance/(sina*sina - cost*cost);
791         pt1.SetXYZ(center.XYZ() - deltacenter*axex);
792         pt2 = pt1;
793         dir1 = normp;
794         dir2.SetXYZ(axex);
795         param1    = param2 = cost*sina*cosa*distance /(sina*sina-cost*cost);
796         param1bis = param2bis = cost*sina*distance / Sqrt(sina*sina-cost*cost);
797
798       }
799       else {                    // on a alors cost > sina
800         typeres = IntAna_Ellipse;
801         nbint = 1;
802         Standard_Real radius = cost*sina*cosa*distance/(cost*cost-sina*sina);
803         deltacenter = sint*sina*sina*distance/(cost*cost-sina*sina);
804         axex.Normalize();
805         pt1.SetXYZ(center.XYZ() + deltacenter*axex);
806         dir1 = normp;
807         dir2.SetXYZ(axex);
808         param1    = radius;
809         param1bis = cost*sina*distance/ Sqrt(cost*cost - sina*sina);
810       }
811     }
812   }
813   
814   //-- On a du mal a gerer plus loin (Value ProjLib, Params ... )
815   //-- des hyperboles trop bizarres
816   //-- On retourne False -> Traitement par biparametree
817   static Standard_Real EllipseLimit   = 1.0E+9; //OCC513(apo) 1000000
818   static Standard_Real HyperbolaLimit = 2.0E+6; //OCC537(apo) 50000
819   if(typeres==IntAna_Ellipse && nbint>=1) { 
820     if(Abs(param1) > EllipseLimit || Abs(param1bis) > EllipseLimit)  { 
821       done=Standard_False; 
822       return;
823     }
824   }
825   if(typeres==IntAna_Hyperbola && nbint>=2) { 
826     if(Abs(param2) > HyperbolaLimit || Abs(param2bis) > HyperbolaLimit)  { 
827       done = Standard_False; 
828       return;
829     }
830   }
831   if(typeres==IntAna_Hyperbola && nbint>=1) { 
832     if(Abs(param1) > HyperbolaLimit || Abs(param1bis) > HyperbolaLimit)  {
833       done=Standard_False;
834       return;
835     }
836   }
837
838   done = Standard_True;
839 }
840
841 //=======================================================================
842 //function : IntAna_QuadQuadGeo
843 //purpose  : Pln Sphere
844 //=======================================================================
845 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P,
846                                        const gp_Sphere& S)
847 : done(Standard_False),
848   nbint(0),
849   typeres(IntAna_Empty),
850   pt1(0,0,0),
851   pt2(0,0,0),
852   pt3(0,0,0),
853   pt4(0,0,0),
854   param1(0),
855   param2(0),
856   param3(0),
857   param4(0),
858   param1bis(0),
859   param2bis(0),
860   myCommonGen(Standard_False),
861   myPChar(0,0,0)
862 {
863   InitTolerances();
864   Perform(P,S);
865 }
866 //=======================================================================
867 //function : Perform
868 //purpose  : 
869 //=======================================================================
870 void IntAna_QuadQuadGeo::Perform( const gp_Pln& P
871                                  ,const gp_Sphere& S) 
872 {
873   
874   done = Standard_False;
875   Standard_Real A,B,C,D,dist, radius;
876   Standard_Real X,Y,Z;
877
878   nbint = 0;
879 // debug JAG : on met typeres = IntAna_Empty par defaut...
880   typeres = IntAna_Empty;
881   
882   P.Coefficients(A,B,C,D);
883   S.Location().Coord(X,Y,Z);
884   radius = S.Radius();
885   
886   dist = A * X + B * Y + C * Z + D;
887   
888   if (Abs( Abs(dist) - radius) < Epsilon(radius)) {
889     // on a une seule solution : le point projection du centre de la sphere
890     // sur le plan
891     nbint = 1;
892     typeres = IntAna_Point;
893     pt1.SetCoord(X - dist*A, Y - dist*B, Z - dist*C);
894   }
895   else if (Abs(dist) < radius) {
896     // on a un cercle solution
897     nbint = 1;
898     typeres = IntAna_Circle;
899     pt1.SetCoord(X - dist*A, Y - dist*B, Z - dist*C);
900     dir1 = P.Axis().Direction();
901     if(P.Direct()==Standard_False) dir1.Reverse();
902     dir2 = P.Position().XDirection();
903     param1 = Sqrt(radius*radius - dist*dist);
904   }
905   param2bis=0.0; //-- pour eviter param2bis not used .... 
906   done = Standard_True;
907 }
908
909 //=======================================================================
910 //function : IntAna_QuadQuadGeo
911 //purpose  : Cylinder - Cylinder
912 //=======================================================================
913 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl1,
914                                        const gp_Cylinder& Cyl2,
915                                        const Standard_Real Tol) 
916 : done(Standard_False),
917   nbint(0),
918   typeres(IntAna_Empty),
919   pt1(0,0,0),
920   pt2(0,0,0),
921   pt3(0,0,0),
922   pt4(0,0,0),
923   param1(0),
924   param2(0),
925   param3(0),
926   param4(0),
927   param1bis(0),
928   param2bis(0),
929   myCommonGen(Standard_False),
930   myPChar(0,0,0)
931 {
932   InitTolerances();
933   Perform(Cyl1,Cyl2,Tol);
934 }
935 //=======================================================================
936 //function : Perform
937 //purpose  : 
938 //=======================================================================
939 void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl1,
940                                  const gp_Cylinder& Cyl2,
941                                  const Standard_Real Tol) 
942 {
943   done=Standard_True;
944   //---------------------------- Parallel axes -------------------------
945   AxeOperator A1A2(Cyl1.Axis(),Cyl2.Axis());
946   Standard_Real R1=Cyl1.Radius();
947   Standard_Real R2=Cyl2.Radius();
948   Standard_Real RmR, RmR_Relative;
949   RmR=(R1>R2)? (R1-R2) : (R2-R1);
950   {
951     Standard_Real Rmax;
952     Rmax=(R1>R2)? R1 : R2;
953     RmR_Relative=RmR/Rmax;
954   }
955
956   Standard_Real DistA1A2=A1A2.Distance();
957   
958   if(A1A2.Parallel())
959   {
960     if(DistA1A2<=Tol)
961     {
962       if(RmR<=Tol)
963       {
964         typeres=IntAna_Same;
965       }
966       else
967       {
968         typeres=IntAna_Empty;
969       }
970     }
971     else
972     {  //-- DistA1A2 > Tol
973       gp_Pnt P1=Cyl1.Location();
974       gp_Pnt P2t=Cyl2.Location();
975       gp_Pnt P2;
976       //-- P2t is projected on the plane (P1,DirCylX,DirCylY)
977       gp_Dir DirCyl = Cyl1.Position().Direction();
978       Standard_Real ProjP2OnDirCyl1=gp_Vec(DirCyl).Dot(gp_Vec(P1,P2t));
979       
980       P2.SetCoord(P2t.X() - ProjP2OnDirCyl1*DirCyl.X(),
981                   P2t.Y() - ProjP2OnDirCyl1*DirCyl.Y(),
982                   P2t.Z() - ProjP2OnDirCyl1*DirCyl.Z());
983       //-- 
984       Standard_Real R1pR2=R1+R2;
985       if(DistA1A2>(R1pR2+Tol))
986       { 
987         typeres=IntAna_Empty;
988         nbint=0;
989       }
990       else if(DistA1A2>(R1pR2))
991       {
992         //-- 1 Tangent line -------------------------------------OK
993         typeres=IntAna_Line;
994
995         nbint=1;
996         dir1=DirCyl;
997         Standard_Real R1_R1pR2=R1/R1pR2;
998         pt1.SetCoord( P1.X() + R1_R1pR2 * (P2.X()-P1.X()),
999                       P1.Y() + R1_R1pR2 * (P2.Y()-P1.Y()),
1000                       P1.Z() + R1_R1pR2 * (P2.Z()-P1.Z()));
1001       }
1002       else if(DistA1A2>RmR)
1003       {
1004         //-- 2 lines ---------------------------------------------OK
1005         typeres=IntAna_Line;
1006         nbint=2;
1007         dir1=DirCyl;
1008         gp_Vec P1P2(P1,P2);
1009         gp_Dir DirA1A2=gp_Dir(P1P2);
1010         gp_Dir Ortho_dir1_P1P2 = dir1.Crossed(DirA1A2);
1011         dir2=dir1;
1012         Standard_Real Alpha=0.5*(R1*R1-R2*R2+DistA1A2*DistA1A2)/(DistA1A2);       
1013
1014         //Standard_Real Beta = Sqrt(R1*R1-Alpha*Alpha);
1015         Standard_Real anSqrtArg = R1*R1-Alpha*Alpha;
1016         Standard_Real Beta = (anSqrtArg > 0.) ? Sqrt(anSqrtArg) : 0.;
1017         
1018         if((Beta+Beta)<Tol)
1019         {
1020           nbint=1;
1021           pt1.SetCoord( P1.X() + Alpha*DirA1A2.X()
1022                        ,P1.Y() + Alpha*DirA1A2.Y()
1023                        ,P1.Z() + Alpha*DirA1A2.Z());
1024         }
1025         else
1026         { 
1027           pt1.SetCoord( P1.X() + Alpha*DirA1A2.X() + Beta*Ortho_dir1_P1P2.X(),
1028                         P1.Y() + Alpha*DirA1A2.Y() + Beta*Ortho_dir1_P1P2.Y(),
1029                         P1.Z() + Alpha*DirA1A2.Z() + Beta*Ortho_dir1_P1P2.Z());
1030           
1031           pt2.SetCoord( P1.X() + Alpha*DirA1A2.X() - Beta*Ortho_dir1_P1P2.X(),
1032                         P1.Y() + Alpha*DirA1A2.Y() - Beta*Ortho_dir1_P1P2.Y(),
1033                         P1.Z() + Alpha*DirA1A2.Z() - Beta*Ortho_dir1_P1P2.Z());
1034         }
1035       }
1036       else if(DistA1A2>(RmR-Tol))
1037       {
1038         //-- 1 Tangent ------------------------------------------OK
1039         typeres=IntAna_Line;
1040         nbint=1;
1041         dir1=DirCyl;
1042         Standard_Real R1_RmR=R1/RmR;
1043
1044         if(R1 < R2)
1045           R1_RmR = -R1_RmR;
1046
1047         pt1.SetCoord( P1.X() + R1_RmR * (P2.X()-P1.X()),
1048                       P1.Y() + R1_RmR * (P2.Y()-P1.Y()),
1049                       P1.Z() + R1_RmR * (P2.Z()-P1.Z()));
1050       }
1051       else {
1052         nbint=0;
1053         typeres=IntAna_Empty;
1054       }
1055     }
1056   }
1057   else { //-- No Parallel Axis ---------------------------------OK
1058     if((RmR_Relative<=myEPSILON_CYLINDER_DELTA_RADIUS) 
1059        && (DistA1A2 <= myEPSILON_CYLINDER_DELTA_DISTANCE))
1060     {
1061       //-- PI/2 between the two axis   and   Intersection  
1062       //-- and identical radius
1063       typeres=IntAna_Ellipse;
1064       nbint=2;
1065       gp_Dir DirCyl1=Cyl1.Position().Direction();
1066       gp_Dir DirCyl2=Cyl2.Position().Direction();
1067       pt1=pt2=A1A2.PtIntersect();
1068       
1069       Standard_Real A=DirCyl1.Angle(DirCyl2);
1070       Standard_Real B;
1071       B=Abs(Sin(0.5*(M_PI-A)));
1072       A=Abs(Sin(0.5*A));
1073       
1074       if(A==0.0 || B==0.0)
1075       {
1076         typeres=IntAna_Same;
1077         return;
1078       }
1079       
1080       gp_Vec dircyl1(DirCyl1);gp_Vec dircyl2(DirCyl2);
1081       dir1 = gp_Dir(dircyl1.Added(dircyl2));
1082       dir2 = gp_Dir(dircyl1.Subtracted(dircyl2));
1083         
1084       param2   = Cyl1.Radius() / A;
1085       param1   = Cyl1.Radius() / B;
1086       param2bis= param1bis = Cyl1.Radius();
1087       if(param1 < param1bis)
1088       {
1089         A=param1;
1090         param1=param1bis;
1091         param1bis=A;
1092       }
1093
1094       if(param2 < param2bis)
1095       {
1096         A=param2;
1097         param2=param2bis;
1098         param2bis=A;
1099       }
1100     }
1101     else
1102     {
1103       if(Abs(DistA1A2-Cyl1.Radius()-Cyl2.Radius())<Tol)
1104       {
1105         typeres = IntAna_Point;
1106         Standard_Real d,p1,p2;
1107
1108         gp_Dir D1 = Cyl1.Axis().Direction();
1109         gp_Dir D2 = Cyl2.Axis().Direction();
1110         A1A2.Distance(d,p1,p2);
1111         gp_Pnt P = Cyl1.Axis().Location();
1112         gp_Pnt P1(P.X() - p1*D1.X(),
1113                   P.Y() - p1*D1.Y(),
1114                   P.Z() - p1*D1.Z());
1115         
1116         P = Cyl2.Axis().Location();
1117         
1118         gp_Pnt P2(P.X() - p2*D2.X(),
1119                   P.Y() - p2*D2.Y(),
1120                   P.Z() - p2*D2.Z());
1121         
1122         gp_Vec P1P2(P1,P2);
1123         D1=gp_Dir(P1P2);
1124         p1=Cyl1.Radius();
1125         
1126         pt1.SetCoord(P1.X() + p1*D1.X(),
1127                      P1.Y() + p1*D1.Y(),
1128                      P1.Z() + p1*D1.Z());
1129         nbint = 1;
1130       }
1131       else
1132       {
1133         typeres=IntAna_NoGeometricSolution;
1134       }
1135     }
1136   }
1137 }
1138 //=======================================================================
1139 //function : IntAna_QuadQuadGeo
1140 //purpose  : Cylinder - Cone
1141 //=======================================================================
1142 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
1143                                        const gp_Cone& Con,
1144                                        const Standard_Real Tol) 
1145 : done(Standard_False),
1146   nbint(0),
1147   typeres(IntAna_Empty),
1148   pt1(0,0,0),
1149   pt2(0,0,0),
1150   pt3(0,0,0),
1151   pt4(0,0,0),
1152   param1(0),
1153   param2(0),
1154   param3(0),
1155   param4(0),
1156   param1bis(0),
1157   param2bis(0),
1158   myCommonGen(Standard_False),
1159   myPChar(0,0,0)
1160 {
1161   InitTolerances();
1162   Perform(Cyl,Con,Tol);
1163 }
1164 //=======================================================================
1165 //function : Perform
1166 //purpose  : 
1167 //=======================================================================
1168   void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl,
1169                                    const gp_Cone& Con,
1170                                    const Standard_Real ) 
1171 {
1172   done=Standard_True;
1173   AxeOperator A1A2(Cyl.Axis(),Con.Axis());
1174   if(A1A2.Same()) {
1175     gp_Pnt Pt=Con.Apex();
1176     Standard_Real dist=Cyl.Radius()/(Tan(Con.SemiAngle()));
1177     gp_Dir dir=Cyl.Position().Direction();
1178     pt1.SetCoord( Pt.X() + dist*dir.X()
1179                  ,Pt.Y() + dist*dir.Y()
1180                  ,Pt.Z() + dist*dir.Z());
1181     pt2.SetCoord( Pt.X() - dist*dir.X()
1182                  ,Pt.Y() - dist*dir.Y()
1183                  ,Pt.Z() - dist*dir.Z());
1184     dir1=dir2=dir;
1185     param1=param2=Cyl.Radius();
1186     nbint=2;
1187     typeres=IntAna_Circle;
1188
1189   }
1190   else {
1191     typeres=IntAna_NoGeometricSolution;
1192   }
1193 }
1194 //=======================================================================
1195 //function : 
1196 //purpose  : Cylinder - Sphere
1197 //=======================================================================
1198   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
1199                                          const gp_Sphere& Sph,
1200                                          const Standard_Real Tol) 
1201 : done(Standard_False),
1202   nbint(0),
1203   typeres(IntAna_Empty),
1204   pt1(0,0,0),
1205   pt2(0,0,0),
1206   pt3(0,0,0),
1207   pt4(0,0,0),
1208   param1(0),
1209   param2(0),
1210   param3(0),
1211   param4(0),
1212   param1bis(0),
1213   param2bis(0),
1214   myCommonGen(Standard_False),
1215   myPChar(0,0,0)
1216 {
1217   InitTolerances();
1218   Perform(Cyl,Sph,Tol);
1219 }
1220 //=======================================================================
1221 //function : Perform
1222 //purpose  : 
1223 //=======================================================================
1224   void IntAna_QuadQuadGeo::Perform( const gp_Cylinder& Cyl
1225                                    ,const gp_Sphere& Sph
1226                                    ,const Standard_Real)
1227 {
1228   done=Standard_True;
1229   gp_Pnt Pt=Sph.Location();
1230   AxeOperator A1A2(Cyl.Axis(),Sph.Position().Axis());
1231   if((A1A2.Intersect()  && Pt.Distance(A1A2.PtIntersect())==0.0 )
1232      || (A1A2.Same()))      {
1233     if(Sph.Radius() < Cyl.Radius()) { 
1234       typeres = IntAna_Empty;
1235     }
1236     else { 
1237       Standard_Real dist=Sqrt( Sph.Radius() * Sph.Radius() - Cyl.Radius() * Cyl.Radius() );
1238       gp_Dir dir=Cyl.Position().Direction();
1239       dir1 = dir2 = dir;
1240       typeres=IntAna_Circle;
1241       pt1.SetCoord( Pt.X() + dist*dir.X()
1242                    ,Pt.Y() + dist*dir.Y()
1243                    ,Pt.Z() + dist*dir.Z());
1244       nbint=1;
1245       param1 = Cyl.Radius();
1246       if(dist>RealEpsilon()) {
1247         pt2.SetCoord( Pt.X() - dist*dir.X()
1248                      ,Pt.Y() - dist*dir.Y()
1249                      ,Pt.Z() - dist*dir.Z());
1250         param2=Cyl.Radius();
1251         nbint=2;
1252       }
1253     }
1254   }
1255   else {
1256     typeres=IntAna_NoGeometricSolution; 
1257   }
1258 }
1259
1260 //=======================================================================
1261 //function : IntAna_QuadQuadGeo
1262 //purpose  : Cone - Cone
1263 //=======================================================================
1264   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cone& Con1,
1265                                          const gp_Cone& Con2,
1266                                          const Standard_Real Tol) 
1267 : done(Standard_False),
1268   nbint(0),
1269   typeres(IntAna_Empty),
1270   pt1(0,0,0),
1271   pt2(0,0,0),
1272   pt3(0,0,0),
1273   pt4(0,0,0),
1274   param1(0),
1275   param2(0),
1276   param3(0),
1277   param4(0),
1278   param1bis(0),
1279   param2bis(0),
1280   myCommonGen(Standard_False),
1281   myPChar(0,0,0)
1282 {
1283   InitTolerances();
1284   Perform(Con1,Con2,Tol);
1285 }
1286 //
1287 //=======================================================================
1288 //function : Perform
1289 //purpose  : 
1290 //=======================================================================
1291   void IntAna_QuadQuadGeo::Perform(const gp_Cone& Con1,
1292                                    const gp_Cone& Con2,
1293                                    const Standard_Real Tol) 
1294 {
1295   done=Standard_True;
1296   //
1297   Standard_Real tg1, tg2, aDA1A2, aTol2;
1298   gp_Pnt aPApex1, aPApex2;
1299
1300   Standard_Real TOL_APEX_CONF = 1.e-10;
1301   
1302   //
1303   tg1=Tan(Con1.SemiAngle());
1304   tg2=Tan(Con2.SemiAngle());
1305
1306   if((tg1 * tg2) < 0.) {
1307     tg2 = -tg2;
1308   }
1309   //
1310   aTol2=Tol*Tol;
1311   aPApex1=Con1.Apex();
1312   aPApex2=Con2.Apex();
1313   aDA1A2=aPApex1.SquareDistance(aPApex2);
1314   //
1315   AxeOperator A1A2(Con1.Axis(),Con2.Axis());
1316   //
1317   // 1
1318   if(A1A2.Same()) {
1319     //-- two circles 
1320     Standard_Real x;
1321     gp_Pnt P=Con1.Apex();
1322     gp_Dir D=Con1.Position().Direction();
1323     Standard_Real d=gp_Vec(D).Dot(gp_Vec(P,Con2.Apex()));
1324     
1325     if(Abs(tg1-tg2)>myEPSILON_ANGLE_CONE) { 
1326       if (fabs(d) < TOL_APEX_CONF) {
1327         typeres = IntAna_Point;
1328         nbint = 1;
1329         pt1 = P;
1330         return;
1331       }
1332       x=(d*tg2)/(tg1+tg2);
1333       pt1.SetCoord( P.X() + x*D.X()
1334                    ,P.Y() + x*D.Y()
1335                    ,P.Z() + x*D.Z());
1336       param1=Abs(x*tg1);
1337
1338       x=(d*tg2)/(tg2-tg1);
1339       pt2.SetCoord( P.X() + x*D.X()
1340                    ,P.Y() + x*D.Y()
1341                    ,P.Z() + x*D.Z());
1342       param2=Abs(x*tg1);
1343       dir1 = dir2 = D;
1344       nbint=2;
1345       typeres=IntAna_Circle;
1346     }
1347     else {
1348       if (fabs(d) < TOL_APEX_CONF) { 
1349         typeres=IntAna_Same;
1350       }
1351       else {
1352         typeres=IntAna_Circle;
1353         nbint=1;
1354         x=d*0.5;
1355         pt1.SetCoord( P.X() + x*D.X()
1356                      ,P.Y() + x*D.Y()
1357                      ,P.Z() + x*D.Z());
1358         param1 = Abs(x * tg1);
1359         dir1 = D;
1360       }
1361     }
1362   } //-- fin A1A2.Same
1363   // 2
1364   else if((Abs(tg1-tg2)<myEPSILON_ANGLE_CONE) && (A1A2.Parallel())) {
1365     //-- voir AnVer12mai98
1366     Standard_Real DistA1A2=A1A2.Distance();
1367     gp_Dir DA1=Con1.Position().Direction();
1368     gp_Vec O1O2(Con1.Apex(),Con2.Apex());
1369     gp_Dir O1O2n(O1O2); // normalization of the vector before projection
1370     Standard_Real O1O2_DA1=gp_Vec(DA1).Dot(gp_Vec(O1O2n));
1371
1372     gp_Vec O1_Proj_A2(O1O2n.X()-O1O2_DA1*DA1.X(),
1373                       O1O2n.Y()-O1O2_DA1*DA1.Y(),
1374                       O1O2n.Z()-O1O2_DA1*DA1.Z());
1375     gp_Dir DB1=gp_Dir(O1_Proj_A2);
1376
1377     Standard_Real yO1O2=O1O2.Dot(gp_Vec(DA1));
1378     Standard_Real ABSTG1 = Abs(tg1);
1379     Standard_Real X2 = (DistA1A2/ABSTG1 - yO1O2)*0.5;
1380     Standard_Real X1 = X2+yO1O2;
1381     
1382     gp_Pnt P1(Con1.Apex().X() + X1*( DA1.X() + ABSTG1*DB1.X()),
1383               Con1.Apex().Y() + X1*( DA1.Y() + ABSTG1*DB1.Y()), 
1384               Con1.Apex().Z() + X1*( DA1.Z() + ABSTG1*DB1.Z()));
1385
1386     gp_Pnt MO1O2(0.5*(Con1.Apex().X()+Con2.Apex().X()),
1387                  0.5*(Con1.Apex().Y()+Con2.Apex().Y()),
1388                  0.5*(Con1.Apex().Z()+Con2.Apex().Z()));
1389     gp_Vec P1MO1O2(P1,MO1O2);
1390     
1391     gp_Dir DA1_X_DB1=DA1.Crossed(DB1);
1392     gp_Dir OrthoPln =  DA1_X_DB1.Crossed(gp_Dir(P1MO1O2));
1393     
1394     IntAna_QuadQuadGeo INTER_QUAD_PLN(gp_Pln(P1,OrthoPln),Con1,Tol,Tol);
1395       if(INTER_QUAD_PLN.IsDone()) {
1396       switch(INTER_QUAD_PLN.TypeInter()) {
1397       case IntAna_Ellipse:         {
1398         typeres=IntAna_Ellipse;
1399         gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
1400         pt1 = E.Location();
1401         dir1 = E.Position().Direction();
1402         dir2 = E.Position().XDirection();
1403         param1 = E.MajorRadius();
1404         param1bis = E.MinorRadius();
1405         nbint = 1;
1406         break; 
1407       }
1408       case IntAna_Circle: {
1409         typeres=IntAna_Circle;
1410         gp_Circ C=INTER_QUAD_PLN.Circle(1);
1411         pt1 = C.Location();
1412         dir1 = C.Position().XDirection();
1413         dir2 = C.Position().YDirection();
1414         param1 = C.Radius();
1415         nbint = 1;
1416         break;
1417       }
1418       case IntAna_Hyperbola: {
1419         typeres=IntAna_Hyperbola;
1420         gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
1421         pt1 = pt2 = H.Location();
1422         dir1 = H.Position().Direction();
1423         dir2 = H.Position().XDirection();
1424         param1 = param2 = H.MajorRadius();
1425         param1bis = param2bis = H.MinorRadius();
1426         nbint = 2;
1427         break;
1428       }
1429       case IntAna_Line: {
1430         typeres=IntAna_Line;
1431         gp_Lin H=INTER_QUAD_PLN.Line(1);
1432         pt1 = pt2 = H.Location();
1433         dir1 = dir2 = H.Position().Direction();
1434         param1 = param2 = 0.0;
1435         param1bis = param2bis = 0.0;
1436         nbint = 2;
1437         break;
1438       }
1439       default:
1440         typeres=IntAna_NoGeometricSolution; 
1441       }
1442     }
1443   }// else if((Abs(tg1-tg2)<EPSILON_ANGLE_CONE) && (A1A2.Parallel()))
1444   // 3
1445   else if (aDA1A2<aTol2) {
1446     //
1447     // When apices are coinsided there can be 3 possible cases
1448     // 3.1 - empty solution (iRet=0)
1449     // 3.2 - one line  when cone1 touches cone2 (iRet=1)
1450     // 3.3 - two lines when cone1 intersects cone2 (iRet=2)
1451     //
1452     Standard_Integer iRet;
1453     Standard_Real aGamma, aBeta1, aBeta2;
1454     Standard_Real aD1, aR1, aTgBeta1, aTgBeta2, aHalfPI;
1455     Standard_Real aCosGamma, aSinGamma, aDx, aR2, aRD2, aD2;
1456     gp_Pnt2d aP0, aPA1, aP1, aPA2;
1457     gp_Vec2d aVAx2;
1458     gp_Ax1 aAx1, aAx2;
1459     //
1460     // Preliminary analysis. Determination of iRet
1461     //
1462     iRet=0;
1463     aHalfPI=0.5*M_PI;
1464     aD1=1.;
1465     aPA1.SetCoord(aD1, 0.);
1466     aP0.SetCoord(0., 0.);
1467     //
1468     aAx1=Con1.Axis();
1469     aAx2=Con2.Axis();
1470     aGamma=aAx1.Angle(aAx2);
1471     if (aGamma>aHalfPI){
1472       aGamma=M_PI-aGamma;
1473     }
1474     aCosGamma=Cos(aGamma);
1475     aSinGamma=Sin(aGamma);
1476     //
1477     aBeta1=Con1.SemiAngle();
1478     aTgBeta1=Tan(aBeta1);
1479     aTgBeta1=Abs(aTgBeta1);
1480     //
1481     aBeta2=Con2.SemiAngle();
1482     aTgBeta2=Tan(aBeta2);
1483     aTgBeta2=Abs(aTgBeta2);
1484     //
1485     aR1=aD1*aTgBeta1;
1486     aP1.SetCoord(aD1, aR1);
1487     //
1488     // PA2
1489     aVAx2.SetCoord(aCosGamma, aSinGamma);
1490     gp_Dir2d aDAx2(aVAx2);
1491     gp_Lin2d aLAx2(aP0, aDAx2);
1492     //
1493     gp_Vec2d aV(aP0, aP1);
1494     aDx=aVAx2.Dot(aV);
1495     aPA2=aP0.Translated(aDx*aDAx2);
1496     //
1497     // aR2
1498     aDx=aPA2.Distance(aP0);
1499     aR2=aDx*aTgBeta2;
1500     //
1501     // aRD2
1502     aRD2=aPA2.Distance(aP1);
1503     //
1504     if (aRD2>(aR2+Tol)) {
1505       iRet=0;
1506       typeres=IntAna_Empty; //nothing
1507       return;
1508     }
1509     //
1510     iRet=1; //touch case => 1 line
1511     if (aRD2<(aR2-Tol)) {
1512       iRet=2;//intersection => couple of lines
1513     }
1514     //
1515     // Finding the solution in 3D
1516     //
1517     Standard_Real aDa;
1518     gp_Pnt aQApex1, aQA1, aQA2, aQX, aQX1, aQX2;
1519     gp_Dir aD3Ax1, aD3Ax2;
1520     gp_Lin aLin;
1521     IntAna_QuadQuadGeo aIntr;
1522     //
1523     aQApex1=Con1.Apex();
1524     aD3Ax1=aAx1.Direction(); 
1525     aQA1.SetCoord(aQApex1.X()+aD1*aD3Ax1.X(),
1526                   aQApex1.Y()+aD1*aD3Ax1.Y(),
1527                   aQApex1.Z()+aD1*aD3Ax1.Z());
1528     //
1529     aDx=aD3Ax1.Dot(aAx2.Direction());
1530     if (aDx<0.) {
1531       aAx2.Reverse();
1532     }
1533     aD3Ax2=aAx2.Direction();
1534     //
1535     aD2=aD1*sqrt((1.+aTgBeta1*aTgBeta1)/(1.+aTgBeta2*aTgBeta2));
1536     //
1537     aQA2.SetCoord(aQApex1.X()+aD2*aD3Ax2.X(),
1538                   aQApex1.Y()+aD2*aD3Ax2.Y(),
1539                   aQApex1.Z()+aD2*aD3Ax2.Z());
1540     //
1541     gp_Pln aPln1(aQA1, aD3Ax1);
1542     gp_Pln aPln2(aQA2, aD3Ax2);
1543     //
1544     aIntr.Perform(aPln1, aPln2, Tol, Tol);
1545     if (!aIntr.IsDone() || 0 == aIntr.NbSolutions()) {
1546       iRet=-1; // just in case. it must not be so
1547       typeres=IntAna_NoGeometricSolution; 
1548       return;
1549     }
1550     //
1551     aLin=aIntr.Line(1);
1552     const gp_Dir& aDLin=aLin.Direction();
1553     gp_Vec aVLin(aDLin);
1554     gp_Pnt aOrig=aLin.Location();
1555     gp_Vec aVr(aQA1, aOrig);
1556     aDx=aVLin.Dot(aVr);
1557     aQX=aOrig.Translated(aDx*aVLin);
1558     //
1559     // Final part
1560     //
1561     typeres=IntAna_Line; 
1562     //
1563     param1=0.;
1564     param2 =0.;
1565     param1bis=0.;
1566     param2bis=0.;
1567     //
1568     if (iRet==1) {
1569       // one line
1570       nbint=1;
1571       pt1=aQApex1;
1572       gp_Vec aVX(aQApex1, aQX);
1573       dir1=gp_Dir(aVX);
1574     }
1575     
1576     else {//iRet=2 
1577       // two lines
1578       nbint=2;
1579       aDa=aQA1.Distance(aQX);
1580       aDx=sqrt(aR1*aR1-aDa*aDa);
1581       aQX1=aQX.Translated(aDx*aVLin);
1582       aQX2=aQX.Translated(-aDx*aVLin);
1583       //
1584       pt1=aQApex1;
1585       pt2=aQApex1;
1586       gp_Vec aVX1(aQApex1, aQX1);
1587       dir1=gp_Dir(aVX1);
1588       gp_Vec aVX2(aQApex1, aQX2);
1589       dir2=gp_Dir(aVX2);
1590     }
1591   } //else if (aDA1A2<aTol2) {
1592   //Case when cones have common generatrix
1593   else if(A1A2.Intersect()) {
1594     //Check if apex of one cone belongs another one
1595     Standard_Real u, v, tol2 = Tol*Tol;
1596     ElSLib::Parameters(Con2, aPApex1, u, v);
1597     gp_Pnt p = ElSLib::Value(u, v, Con2);
1598     if(aPApex1.SquareDistance(p) > tol2) {
1599       typeres=IntAna_NoGeometricSolution; 
1600       return;
1601     }
1602     //
1603     ElSLib::Parameters(Con1, aPApex2, u, v);
1604     p = ElSLib::Value(u, v, Con1);
1605     if(aPApex2.SquareDistance(p) > tol2) {
1606       typeres=IntAna_NoGeometricSolution; 
1607       return;
1608     }
1609
1610     //Cones have a common generatrix passing through apexes
1611     myCommonGen = Standard_True;
1612
1613     //common generatrix of cones
1614     gp_Lin aGen(aPApex1, gp_Dir(gp_Vec(aPApex1, aPApex2)));
1615
1616     //Intersection point of axes
1617     gp_Pnt aPAxeInt = A1A2.PtIntersect();
1618
1619     //Characteristic point of intersection curve
1620     u = ElCLib::Parameter(aGen, aPAxeInt);
1621     myPChar = ElCLib::Value(u, aGen);
1622
1623
1624     //Other generatrixes of cones laying in maximal plane
1625     gp_Lin aGen1 = aGen.Rotated(Con1.Axis(), M_PI); 
1626     gp_Lin aGen2 = aGen.Rotated(Con2.Axis(), M_PI); 
1627     //
1628     //Intersection point of generatrixes
1629     gp_Dir aN; //solution plane normal
1630     gp_Dir aD1 = aGen1.Direction();
1631
1632     gp_Dir aD2(aD1.Crossed(aGen.Direction()));
1633
1634     if(aD1.IsParallel(aGen2.Direction(), Precision::Angular())) {
1635       aN = aD1.Crossed(aD2);
1636     }
1637     else if(aGen1.SquareDistance(aGen2) > tol2) {
1638       //Something wrong ???
1639       typeres=IntAna_NoGeometricSolution; 
1640       return;
1641     }
1642     else {
1643       gp_Dir D1 = aGen1.Position().Direction();
1644       gp_Dir D2 = aGen2.Position().Direction();
1645       gp_Pnt O1 = aGen1.Location();
1646       gp_Pnt O2 = aGen2.Location();
1647       Standard_Real D1DotD2 = D1.Dot(D2);
1648       Standard_Real aSin = 1.-D1DotD2*D1DotD2;
1649       gp_Vec O1O2 (O1,O2);
1650       Standard_Real U2 = (D1.XYZ()*(O1O2.Dot(D1))-(O1O2.XYZ())).Dot(D2.XYZ());
1651       U2 /= aSin;
1652       gp_Pnt aPGint(ElCLib::Value(U2, aGen2));
1653     
1654       aD1 = gp_Dir(gp_Vec(aPGint, myPChar));
1655       aN = aD1.Crossed(aD2);
1656     }
1657     //Plane that must contain intersection curves
1658     gp_Pln anIntPln(myPChar, aN);
1659
1660     IntAna_QuadQuadGeo INTER_QUAD_PLN(anIntPln,Con1,Tol,Tol);
1661
1662       if(INTER_QUAD_PLN.IsDone()) {
1663       switch(INTER_QUAD_PLN.TypeInter()) {
1664       case IntAna_Ellipse:         {
1665         typeres=IntAna_Ellipse;
1666         gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
1667         pt1 = E.Location();
1668         dir1 = E.Position().Direction();
1669         dir2 = E.Position().XDirection();
1670         param1 = E.MajorRadius();
1671         param1bis = E.MinorRadius();
1672         nbint = 1;
1673         break; 
1674       }
1675       case IntAna_Circle: {
1676         typeres=IntAna_Circle;
1677         gp_Circ C=INTER_QUAD_PLN.Circle(1);
1678         pt1 = C.Location();
1679         dir1 = C.Position().XDirection();
1680         dir2 = C.Position().YDirection();
1681         param1 = C.Radius();
1682         nbint = 1;
1683         break;
1684       }
1685       case IntAna_Parabola: {
1686         typeres=IntAna_Parabola;
1687         gp_Parab Prb=INTER_QUAD_PLN.Parabola(1);
1688         pt1 = Prb.Location();
1689         dir1 = Prb.Position().Direction();
1690         dir2 = Prb.Position().XDirection();
1691         param1 = Prb.Focal();
1692         nbint = 1;
1693         break;
1694       }
1695       case IntAna_Hyperbola: {
1696         typeres=IntAna_Hyperbola;
1697         gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
1698         pt1 = pt2 = H.Location();
1699         dir1 = H.Position().Direction();
1700         dir2 = H.Position().XDirection();
1701         param1 = param2 = H.MajorRadius();
1702         param1bis = param2bis = H.MinorRadius();
1703         nbint = 2;
1704         break;
1705       }
1706       default:
1707         typeres=IntAna_NoGeometricSolution; 
1708       }
1709     }
1710   }
1711   
1712   else {
1713     typeres=IntAna_NoGeometricSolution; 
1714   }
1715 }
1716 //=======================================================================
1717 //function : IntAna_QuadQuadGeo
1718 //purpose  : Sphere - Cone
1719 //=======================================================================
1720   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Sphere& Sph,
1721                                          const gp_Cone& Con,
1722                                          const Standard_Real Tol) 
1723 : done(Standard_False),
1724   nbint(0),
1725   typeres(IntAna_Empty),
1726   pt1(0,0,0),
1727   pt2(0,0,0),
1728   pt3(0,0,0),
1729   pt4(0,0,0),
1730   param1(0),
1731   param2(0),
1732   param3(0),
1733   param4(0),
1734   param1bis(0),
1735   param2bis(0),
1736   myCommonGen(Standard_False),
1737   myPChar(0,0,0)
1738 {
1739   InitTolerances();
1740   Perform(Sph,Con,Tol);
1741 }
1742 //=======================================================================
1743 //function : Perform
1744 //purpose  : 
1745 //=======================================================================
1746   void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph,
1747                                    const gp_Cone& Con,
1748                                    const Standard_Real)
1749 {
1750   
1751   //
1752   done=Standard_True;
1753   //
1754   AxeOperator A1A2(Con.Axis(),Sph.Position().Axis());
1755   gp_Pnt Pt=Sph.Location();
1756   //
1757   if((A1A2.Intersect() && (Pt.Distance(A1A2.PtIntersect())==0.0))
1758      || A1A2.Same()) {
1759     gp_Pnt ConApex= Con.Apex();
1760     Standard_Real dApexSphCenter=Pt.Distance(ConApex); 
1761     gp_Dir ConDir;
1762     if(dApexSphCenter>RealEpsilon()) { 
1763       ConDir = gp_Dir(gp_Vec(ConApex,Pt));
1764     }
1765     else { 
1766       ConDir = Con.Position().Direction();
1767     }
1768     
1769     Standard_Real Rad=Sph.Radius();
1770     Standard_Real tga=Tan(Con.SemiAngle());
1771
1772
1773     //-- 2 circles
1774     //-- x: Roots of    (x**2 + y**2 = Rad**2)
1775     //--                tga = y / (x+dApexSphCenter)
1776     Standard_Real tgatga = tga * tga;
1777     math_DirectPolynomialRoots Eq( 1.0+tgatga
1778                                   ,2.0*tgatga*dApexSphCenter
1779                                   ,-Rad*Rad + dApexSphCenter*dApexSphCenter*tgatga);
1780     if(Eq.IsDone()) {
1781       Standard_Integer nbsol=Eq.NbSolutions();
1782       if(nbsol==0) {
1783         typeres=IntAna_Empty;
1784       }
1785       else { 
1786         typeres=IntAna_Circle;
1787         if(nbsol>=1) {
1788           Standard_Real x = Eq.Value(1);
1789           Standard_Real dApexSphCenterpx = dApexSphCenter+x;
1790           nbint=1;
1791           pt1.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
1792                        ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
1793                        ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
1794           param1 = tga * dApexSphCenterpx;
1795           param1 = Abs(param1);
1796           dir1 = ConDir;
1797           if(param1<=myEPSILON_MINI_CIRCLE_RADIUS) {
1798             typeres=IntAna_PointAndCircle;
1799             param1=0.0;
1800           }
1801         }
1802         if(nbsol>=2) {
1803           Standard_Real x=Eq.Value(2);
1804           Standard_Real dApexSphCenterpx = dApexSphCenter+x;
1805           nbint=2;
1806           pt2.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
1807                        ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
1808                        ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
1809           param2 = tga * dApexSphCenterpx;
1810           param2 = Abs(param2);
1811           dir2=ConDir;
1812           if(param2<=myEPSILON_MINI_CIRCLE_RADIUS) {
1813             typeres=IntAna_PointAndCircle;
1814             param2=0.0;
1815           }
1816         }
1817       }
1818     }
1819     else {
1820       done=Standard_False;
1821     }
1822   }
1823   else {
1824     typeres=IntAna_NoGeometricSolution; 
1825   }
1826 }
1827
1828 //=======================================================================
1829 //function : IntAna_QuadQuadGeo
1830 //purpose  : Sphere - Sphere
1831 //=======================================================================
1832   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo( const gp_Sphere& Sph1
1833                                          ,const gp_Sphere& Sph2
1834                                          ,const Standard_Real Tol) 
1835 : done(Standard_False),
1836   nbint(0),
1837   typeres(IntAna_Empty),
1838   pt1(0,0,0),
1839   pt2(0,0,0),
1840   pt3(0,0,0),
1841   pt4(0,0,0),
1842   param1(0),
1843   param2(0),
1844   param3(0),
1845   param4(0),
1846   param1bis(0),
1847   param2bis(0),
1848   myCommonGen(Standard_False),
1849   myPChar(0,0,0)
1850 {
1851   InitTolerances();
1852   Perform(Sph1,Sph2,Tol);
1853 }
1854 //=======================================================================
1855 //function : Perform
1856 //purpose  : 
1857 //=======================================================================
1858   void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph1,
1859                                    const gp_Sphere& Sph2,
1860                                    const Standard_Real Tol)   
1861 {
1862   done=Standard_True;
1863   gp_Pnt O1=Sph1.Location();
1864   gp_Pnt O2=Sph2.Location();
1865   Standard_Real dO1O2=O1.Distance(O2);
1866   Standard_Real R1=Sph1.Radius();
1867   Standard_Real R2=Sph2.Radius();
1868   Standard_Real Rmin,Rmax;
1869   typeres=IntAna_Empty;
1870   param2bis=0.0; //-- pour eviter param2bis not used .... 
1871
1872   if(R1>R2) { Rmin=R2; Rmax=R1; } else { Rmin=R1; Rmax=R2; }
1873   
1874   if(dO1O2<=Tol && (Abs(R1-R2) <= Tol)) {
1875     typeres = IntAna_Same;
1876   }
1877   else { 
1878     if(dO1O2<=Tol) { return; } 
1879     gp_Dir Dir=gp_Dir(gp_Vec(O1,O2));
1880     Standard_Real t = Rmax - dO1O2 - Rmin;
1881
1882     //----------------------------------------------------------------------
1883     //--        |----------------- R1 --------------------|
1884     //--        |----dO1O2-----|-----------R2----------|
1885     //--                                            --->--<-- t
1886     //--
1887     //--        |------ R1 ------|---------dO1O2----------|
1888     //--     |-------------------R2-----------------------|
1889     //--  --->--<-- t
1890     //----------------------------------------------------------------------
1891     if(t >= 0.0  && t <=Tol) { 
1892       typeres = IntAna_Point;
1893       nbint = 1;
1894       Standard_Real t2;
1895       if(R1==Rmax) t2=(R1 + (R2 + dO1O2)) * 0.5;
1896       else         t2=(-R1+(dO1O2-R2))*0.5;
1897         
1898       pt1.SetCoord( O1.X() + t2*Dir.X()
1899                    ,O1.Y() + t2*Dir.Y()
1900                    ,O1.Z() + t2*Dir.Z());
1901     }
1902     else  {
1903       //-----------------------------------------------------------------
1904       //--        |----------------- dO1O2 --------------------|
1905       //--        |----R1-----|-----------R2----------|-Tol-|
1906       //--                                            
1907       //--        |----------------- Rmax --------------------|
1908       //--        |----Rmin----|-------dO1O2-------|-Tol-|
1909       //--                                            
1910       //-----------------------------------------------------------------
1911       if((dO1O2 > (R1+R2+Tol)) || (Rmax > (dO1O2+Rmin+Tol))) {
1912         typeres=IntAna_Empty;
1913       }
1914       else {
1915         //---------------------------------------------------------------
1916         //--     
1917         //--
1918         //---------------------------------------------------------------
1919         Standard_Real Alpha=0.5*(R1*R1-R2*R2+dO1O2*dO1O2)/(dO1O2);       
1920         Standard_Real Beta = R1*R1-Alpha*Alpha;
1921         Beta = (Beta>0.0)? Sqrt(Beta) : 0.0;
1922         
1923         if(Beta<= myEPSILON_MINI_CIRCLE_RADIUS) { 
1924           typeres = IntAna_Point;
1925           Alpha = (R1 + (dO1O2 - R2)) * 0.5;
1926         }
1927         else { 
1928           typeres = IntAna_Circle;
1929           dir1 = Dir;
1930           param1 = Beta;
1931         }          
1932         pt1.SetCoord( O1.X() + Alpha*Dir.X()
1933                      ,O1.Y() + Alpha*Dir.Y()
1934                      ,O1.Z() + Alpha*Dir.Z());
1935         
1936         nbint=1;
1937       }
1938     }
1939   }
1940 }
1941
1942 //=======================================================================
1943 //function : IntAna_QuadQuadGeo
1944 //purpose  : Plane - Torus
1945 //=======================================================================
1946 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& Pln,
1947                                        const gp_Torus& Tor,
1948                                        const Standard_Real Tol) 
1949 : done(Standard_False),
1950   nbint(0),
1951   typeres(IntAna_Empty),
1952   pt1(0,0,0),
1953   pt2(0,0,0),
1954   pt3(0,0,0),
1955   pt4(0,0,0),
1956   param1(0),
1957   param2(0),
1958   param3(0),
1959   param4(0),
1960   param1bis(0),
1961   param2bis(0),
1962   myCommonGen(Standard_False),
1963   myPChar(0,0,0)
1964 {
1965   InitTolerances();
1966   Perform(Pln,Tor,Tol);
1967 }
1968 //=======================================================================
1969 //function : Perform
1970 //purpose  : 
1971 //=======================================================================
1972 void IntAna_QuadQuadGeo::Perform(const gp_Pln& Pln,
1973                                  const gp_Torus& Tor,
1974                                  const Standard_Real Tol)
1975 {
1976   done = Standard_True;
1977   //
1978   Standard_Real aRMin, aRMaj;
1979   //
1980   aRMin = Tor.MinorRadius();
1981   aRMaj = Tor.MajorRadius();
1982   if (aRMin >= aRMaj) {
1983     typeres = IntAna_NoGeometricSolution; 
1984     return;
1985   }
1986   //
1987   const gp_Ax1 aPlnAx = Pln.Axis();
1988   const gp_Ax1 aTorAx = Tor.Axis();
1989   //
1990   Standard_Boolean bParallel, bNormal;
1991   //
1992   bParallel = aTorAx.IsParallel(aPlnAx, myEPSILON_AXES_PARA);
1993   bNormal = !bParallel ? aTorAx.IsNormal(aPlnAx, myEPSILON_AXES_PARA) : Standard_False;
1994   if (!bNormal && !bParallel) {
1995     typeres = IntAna_NoGeometricSolution; 
1996     return;
1997   }
1998   //
1999   Standard_Real aDist;
2000   //
2001   gp_Pnt aTorLoc = aTorAx.Location();
2002   if (bParallel) {
2003     Standard_Real aDt, X, Y, Z, A, B, C, D;
2004     //
2005     Pln.Coefficients(A,B,C,D);
2006     aTorLoc.Coord(X,Y,Z);
2007     aDist = A*X + B*Y + C*Z + D;
2008     //
2009     if ((Abs(aDist) - aRMin) > Tol) {
2010       typeres=IntAna_Empty;
2011       return;
2012     }
2013     //
2014     typeres = IntAna_Circle;
2015     //
2016     pt1.SetCoord(X - aDist*A, Y - aDist*B, Z - aDist*C);
2017     aDt = Sqrt(Abs(aRMin*aRMin - aDist*aDist));
2018     param1 = aRMaj + aDt;
2019     dir1 = aTorAx.Direction();
2020     nbint = 1;
2021     if ((Abs(aDist) < aRMin) && (aDt > Tol)) {
2022       pt2 = pt1;
2023       param2 = aRMaj - aDt;
2024       dir2 = dir1;
2025       nbint = 2;
2026     }
2027   }
2028   //
2029   else {
2030     aDist = Pln.Distance(aTorLoc);
2031     if (aDist > myEPSILON_DISTANCE) {
2032       typeres = IntAna_NoGeometricSolution; 
2033       return;
2034     }
2035     //
2036     typeres = IntAna_Circle;
2037     param2 = param1 = aRMin;
2038     dir2 = dir1 = aPlnAx.Direction();
2039     nbint = 2;
2040     //
2041     gp_Dir aDir = aTorAx.Direction()^dir1;
2042     pt1.SetXYZ(aTorLoc.XYZ() + aRMaj*aDir.XYZ());
2043     pt2.SetXYZ(aTorLoc.XYZ() - aRMaj*aDir.XYZ());
2044   }
2045 }
2046
2047 //=======================================================================
2048 //function : IntAna_QuadQuadGeo
2049 //purpose  : Cylinder - Torus
2050 //=======================================================================
2051 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
2052                                        const gp_Torus& Tor,
2053                                        const Standard_Real Tol) 
2054 : done(Standard_False),
2055   nbint(0),
2056   typeres(IntAna_Empty),
2057   pt1(0,0,0),
2058   pt2(0,0,0),
2059   pt3(0,0,0),
2060   pt4(0,0,0),
2061   param1(0),
2062   param2(0),
2063   param3(0),
2064   param4(0),
2065   param1bis(0),
2066   param2bis(0),
2067   myCommonGen(Standard_False),
2068   myPChar(0,0,0)
2069 {
2070   InitTolerances();
2071   Perform(Cyl,Tor,Tol);
2072 }
2073 //=======================================================================
2074 //function : Perform
2075 //purpose  : 
2076 //=======================================================================
2077 void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl,
2078                                  const gp_Torus& Tor,
2079                                  const Standard_Real Tol) 
2080 {
2081   done = Standard_True;
2082   //
2083   Standard_Real aRMin, aRMaj;
2084   //
2085   aRMin = Tor.MinorRadius();
2086   aRMaj = Tor.MajorRadius();
2087   if (aRMin >= aRMaj) {
2088     typeres = IntAna_NoGeometricSolution; 
2089     return;
2090   }
2091   //
2092   const gp_Ax1 aCylAx = Cyl.Axis();
2093   const gp_Ax1 aTorAx = Tor.Axis();
2094   //
2095   const gp_Lin aLin(aTorAx);
2096   const gp_Pnt aLocCyl = Cyl.Location();
2097   //
2098   if (!aTorAx.IsParallel(aCylAx, myEPSILON_AXES_PARA) ||
2099       (aLin.Distance(aLocCyl) > myEPSILON_DISTANCE)) {
2100     typeres = IntAna_NoGeometricSolution; 
2101     return;
2102   }
2103   //
2104   Standard_Real aRCyl;
2105   //
2106   aRCyl = Cyl.Radius();
2107   if (((aRCyl + Tol) < (aRMaj - aRMin)) || ((aRCyl - Tol) > (aRMaj + aRMin))) {
2108     typeres = IntAna_Empty;
2109     return;
2110   }
2111   //
2112   typeres = IntAna_Circle;
2113   //
2114   Standard_Real aDist = Sqrt(Abs(aRMin*aRMin - (aRCyl-aRMaj)*(aRCyl-aRMaj)));
2115   gp_XYZ aTorLoc = aTorAx.Location().XYZ();
2116   //
2117   dir1 = aTorAx.Direction();
2118   pt1.SetXYZ(aTorLoc + aDist*dir1.XYZ());
2119   param1 = aRCyl;
2120   nbint = 1;
2121   if ((aDist > Tol) && (aRCyl > (aRMaj - aRMin)) &&
2122                        (aRCyl < (aRMaj + aRMin))) {
2123     dir2 = dir1;
2124     pt2.SetXYZ(aTorLoc - aDist*dir2.XYZ());
2125     param2 = param1;
2126     nbint = 2;
2127   }
2128 }
2129
2130 //=======================================================================
2131 //function : IntAna_QuadQuadGeo
2132 //purpose  : Cone - Torus
2133 //=======================================================================
2134 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cone& Con,
2135                                        const gp_Torus& Tor,
2136                                        const Standard_Real Tol) 
2137 : done(Standard_False),
2138   nbint(0),
2139   typeres(IntAna_Empty),
2140   pt1(0,0,0),
2141   pt2(0,0,0),
2142   pt3(0,0,0),
2143   pt4(0,0,0),
2144   param1(0),
2145   param2(0),
2146   param3(0),
2147   param4(0),
2148   param1bis(0),
2149   param2bis(0),
2150   myCommonGen(Standard_False),
2151   myPChar(0,0,0)
2152 {
2153   InitTolerances();
2154   Perform(Con,Tor,Tol);
2155 }
2156 //=======================================================================
2157 //function : Perform
2158 //purpose  : 
2159 //=======================================================================
2160 void IntAna_QuadQuadGeo::Perform(const gp_Cone& Con,
2161                                  const gp_Torus& Tor,
2162                                  const Standard_Real Tol) 
2163 {
2164   done = Standard_True;
2165   //
2166   Standard_Real aRMin, aRMaj;
2167   //
2168   aRMin = Tor.MinorRadius();
2169   aRMaj = Tor.MajorRadius();
2170   if (aRMin >= aRMaj) {
2171     typeres = IntAna_NoGeometricSolution; 
2172     return;
2173   }
2174   //
2175   const gp_Ax1 aConAx = Con.Axis();
2176   const gp_Ax1 aTorAx = Tor.Axis();
2177   //
2178   const gp_Lin aLin(aTorAx);
2179   const gp_Pnt aConApex = Con.Apex();
2180   //
2181   if (!aTorAx.IsParallel(aConAx, myEPSILON_AXES_PARA) ||
2182       (aLin.Distance(aConApex) > myEPSILON_DISTANCE)) {
2183     typeres = IntAna_NoGeometricSolution; 
2184     return;
2185   }
2186   //
2187   Standard_Real anAngle, aDist, aParam[4], aDt;
2188   Standard_Integer i;
2189   gp_Pnt aTorLoc, aPCT, aPN, aPt[4];
2190   gp_Dir aDir[4];
2191   //
2192   anAngle = Con.SemiAngle();
2193   aTorLoc = aTorAx.Location();
2194   //
2195   aPN.SetXYZ(aTorLoc.XYZ() + aRMaj*Tor.YAxis().Direction().XYZ());
2196   gp_Dir aDN (gp_Vec(aTorLoc, aPN));
2197   gp_Ax1 anAxCLRot(aConApex, aDN);
2198   gp_Lin aConL = aLin.Rotated(anAxCLRot, anAngle);
2199   gp_Dir aDL = aConL.Position().Direction();
2200   gp_Dir aXDir = Tor.XAxis().Direction();
2201   //
2202   typeres = IntAna_Empty;
2203   //
2204   for (i = 0; i < 2; ++i) {
2205     if (i) {
2206       aXDir.Reverse();
2207     }
2208     aPCT.SetXYZ(aTorLoc.XYZ() + aRMaj*aXDir.XYZ());
2209     //
2210     aDist = aConL.Distance(aPCT);
2211     if (aDist > aRMin+Tol) {
2212       continue;
2213     }
2214     //
2215     typeres = IntAna_Circle;
2216     //
2217     gp_XYZ aPh = aPCT.XYZ() - aDist*aConL.Normal(aPCT).Direction().XYZ();
2218     aDt = Sqrt(Abs(aRMin*aRMin - aDist*aDist));
2219     //
2220     gp_Pnt aP;
2221     gp_XYZ aDVal = aDt*aDL.XYZ();
2222     aP.SetXYZ(aPh + aDVal);
2223     aParam[nbint] = aLin.Distance(aP);
2224     aPt[nbint].SetXYZ(aP.XYZ() - aParam[nbint]*aXDir.XYZ());
2225     aDir[nbint] = aTorAx.Direction();
2226     ++nbint;
2227     if ((aDist < aRMin) && (aDt > Tol)) {
2228       aP.SetXYZ(aPh - aDVal);
2229       aParam[nbint] = aLin.Distance(aP);
2230       aPt[nbint].SetXYZ(aP.XYZ() - aParam[nbint]*aXDir.XYZ());
2231       aDir[nbint] = aDir[nbint-1];
2232       ++nbint;
2233     }
2234   }
2235   //
2236   for (i = 0; i < nbint; ++i) {
2237     switch (i) {
2238     case 0:{
2239       pt1 = aPt[i];
2240       param1 = aParam[i];
2241       dir1 = aDir[i];
2242       break;
2243     }
2244     case 1:{
2245       pt2 = aPt[i];
2246       param2 = aParam[i];
2247       dir2 = aDir[i];
2248       break;
2249     }
2250     case 2:{
2251       pt3 = aPt[i];
2252       param3 = aParam[i];
2253       dir3 = aDir[i];
2254       break;
2255     }
2256     case 3:{
2257       pt4 = aPt[i];
2258       param4 = aParam[i];
2259       dir4 = aDir[i];
2260       break;
2261     }
2262     default:
2263       break;
2264     }
2265   }
2266 }
2267
2268 //=======================================================================
2269 //function : IntAna_QuadQuadGeo
2270 //purpose  : Sphere - Torus
2271 //=======================================================================
2272 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Sphere& Sph,
2273                                        const gp_Torus& Tor,
2274                                        const Standard_Real Tol) 
2275 : done(Standard_False),
2276   nbint(0),
2277   typeres(IntAna_Empty),
2278   pt1(0,0,0),
2279   pt2(0,0,0),
2280   pt3(0,0,0),
2281   pt4(0,0,0),
2282   param1(0),
2283   param2(0),
2284   param3(0),
2285   param4(0),
2286   param1bis(0),
2287   param2bis(0),
2288   myCommonGen(Standard_False),
2289   myPChar(0,0,0)
2290 {
2291   InitTolerances();
2292   Perform(Sph,Tor,Tol);
2293 }
2294 //=======================================================================
2295 //function : Perform
2296 //purpose  : 
2297 //=======================================================================
2298 void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph,
2299                                  const gp_Torus& Tor,
2300                                  const Standard_Real Tol) 
2301 {
2302   done = Standard_True;
2303   //
2304   Standard_Real aRMin, aRMaj;
2305   //
2306   aRMin = Tor.MinorRadius();
2307   aRMaj = Tor.MajorRadius();
2308   if (aRMin >= aRMaj) {
2309     typeres = IntAna_NoGeometricSolution; 
2310     return;
2311   }
2312   //
2313   const gp_Ax1 aTorAx = Tor.Axis();
2314   const gp_Lin aLin(aTorAx);
2315   const gp_Pnt aSphLoc = Sph.Location();
2316   //
2317   if (aLin.Distance(aSphLoc) > myEPSILON_DISTANCE) {
2318     typeres = IntAna_NoGeometricSolution;
2319     return;
2320   }
2321   //
2322   Standard_Real aRSph, aDist;
2323   gp_Pnt aTorLoc;
2324   //
2325   gp_Dir aXDir = Tor.XAxis().Direction();
2326   aTorLoc.SetXYZ(aTorAx.Location().XYZ() + aRMaj*aXDir.XYZ());
2327   aRSph = Sph.Radius();
2328   //
2329   gp_Vec aVec12(aTorLoc, aSphLoc);
2330   aDist = aVec12.Magnitude();
2331   if (((aDist - Tol) > (aRMin + aRSph)) || 
2332       ((aDist + Tol) < Abs(aRMin - aRSph))) {
2333     typeres = IntAna_Empty;
2334     return;
2335   }
2336   //
2337   typeres = IntAna_Circle;
2338   //
2339   Standard_Real anAlpha, aBeta;
2340   //
2341   anAlpha = 0.5*(aRMin*aRMin - aRSph*aRSph + aDist*aDist ) / aDist;
2342   aBeta = Sqrt(Abs(aRMin*aRMin - anAlpha*anAlpha));
2343   //
2344   gp_Dir aDir12(aVec12);
2345   gp_XYZ aPh = aTorLoc.XYZ() + anAlpha*aDir12.XYZ();
2346   gp_Dir aDC = Tor.YAxis().Direction()^aDir12;
2347   //
2348   gp_Pnt aP;
2349   gp_XYZ aDVal = aBeta*aDC.XYZ();
2350   aP.SetXYZ(aPh + aDVal);
2351   param1 = aLin.Distance(aP);
2352   pt1.SetXYZ(aP.XYZ() - param1*aXDir.XYZ());
2353   dir1 = aTorAx.Direction();
2354   nbint = 1;
2355   if ((aDist < (aRSph + aRMin)) && (aDist > Abs(aRSph - aRMin)) && 
2356       (aDVal.Modulus() > Tol)) {
2357     aP.SetXYZ(aPh - aDVal);
2358     param2 = aLin.Distance(aP);
2359     pt2.SetXYZ(aP.XYZ() - param2*aXDir.XYZ());
2360     dir2 = dir1;
2361     nbint = 2;
2362   }
2363 }
2364
2365 //=======================================================================
2366 //function : IntAna_QuadQuadGeo
2367 //purpose  : Torus - Torus
2368 //=======================================================================
2369 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Torus& Tor1,
2370                                        const gp_Torus& Tor2,
2371                                        const Standard_Real Tol) 
2372 : done(Standard_False),
2373   nbint(0),
2374   typeres(IntAna_Empty),
2375   pt1(0,0,0),
2376   pt2(0,0,0),
2377   pt3(0,0,0),
2378   pt4(0,0,0),
2379   param1(0),
2380   param2(0),
2381   param3(0),
2382   param4(0),
2383   param1bis(0),
2384   param2bis(0),
2385   myCommonGen(Standard_False),
2386   myPChar(0,0,0)
2387 {
2388   InitTolerances();
2389   Perform(Tor1,Tor2,Tol);
2390 }
2391 //=======================================================================
2392 //function : Perform
2393 //purpose  : 
2394 //=======================================================================
2395 void IntAna_QuadQuadGeo::Perform(const gp_Torus& Tor1,
2396                                  const gp_Torus& Tor2,
2397                                  const Standard_Real Tol) 
2398 {
2399   done = Standard_True;
2400   //
2401   Standard_Real aRMin1, aRMin2, aRMaj1, aRMaj2;
2402   //
2403   aRMin1 = Tor1.MinorRadius();
2404   aRMaj1 = Tor1.MajorRadius();
2405   aRMin2 = Tor2.MinorRadius();
2406   aRMaj2 = Tor2.MajorRadius();
2407   if (aRMin1 >= aRMaj1 || aRMin2 >= aRMaj2) {
2408     typeres = IntAna_NoGeometricSolution;
2409     return;
2410   }
2411   //
2412   const gp_Ax1 anAx1 = Tor1.Axis();
2413   const gp_Ax1 anAx2 = Tor2.Axis();
2414   //
2415   gp_Lin aL1(anAx1);
2416   if (!anAx1.IsParallel(anAx2, myEPSILON_AXES_PARA) ||
2417       (aL1.Distance(anAx2.Location()) > myEPSILON_DISTANCE)) {
2418     typeres = IntAna_NoGeometricSolution;
2419     return;
2420   }
2421   //
2422   gp_Pnt aLoc1, aLoc2;
2423   //
2424   aLoc1 = anAx1.Location();
2425   aLoc2 = anAx2.Location();
2426   //
2427   if (aLoc1.IsEqual(aLoc2, Tol) &&
2428       (Abs(aRMin1 - aRMin2) <= Tol) && 
2429       (Abs(aRMaj1 - aRMaj2) <= Tol)) {
2430     typeres = IntAna_Same;
2431     return;
2432   }
2433   //
2434   Standard_Real aDist;
2435   gp_Pnt aP1, aP2;
2436   //
2437   gp_Dir aXDir1 = Tor1.XAxis().Direction();
2438   aP1.SetXYZ(aLoc1.XYZ() + aRMaj1*aXDir1.XYZ());
2439   aP2.SetXYZ(aLoc2.XYZ() + aRMaj2*aXDir1.XYZ());
2440   //
2441   gp_Vec aV12(aP1, aP2);
2442   aDist = aV12.Magnitude();
2443   if (((aDist - Tol) > (aRMin1 + aRMin2)) || 
2444       ((aDist + Tol) < Abs(aRMin1 - aRMin2))) {
2445     typeres = IntAna_Empty;
2446     return;
2447   }
2448   //
2449   typeres = IntAna_Circle;
2450   //
2451   Standard_Real anAlpha, aBeta;
2452   //
2453   anAlpha = 0.5*(aRMin1*aRMin1 - aRMin2*aRMin2 + aDist*aDist ) / aDist;
2454   aBeta = Sqrt(Abs(aRMin1*aRMin1 - anAlpha*anAlpha));
2455   //
2456   gp_Dir aDir12(aV12);
2457   gp_XYZ aPh = aP1.XYZ() + anAlpha*aDir12.XYZ();
2458   gp_Dir aDC = Tor1.YAxis().Direction()^aDir12;
2459   //
2460   gp_Pnt aP;
2461   gp_XYZ aDVal = aBeta*aDC.XYZ();
2462   aP.SetXYZ(aPh + aDVal);
2463   param1 = aL1.Distance(aP);
2464   pt1.SetXYZ(aP.XYZ() - param1*aXDir1.XYZ());
2465   dir1 = anAx1.Direction();
2466   nbint = 1;
2467   if ((aDist < (aRMin1 + aRMin2)) && (aDist > Abs(aRMin1 - aRMin2)) && 
2468       aDVal.Modulus() > Tol) {
2469     aP.SetXYZ(aPh - aDVal);
2470     param2 = aL1.Distance(aP);
2471     pt2.SetXYZ(aP.XYZ() - param2*aXDir1.XYZ());
2472     dir2 = dir1;
2473     nbint = 2;
2474   }
2475 }
2476
2477 //=======================================================================
2478 //function : Point
2479 //purpose  : Returns a Point
2480 //=======================================================================
2481   gp_Pnt IntAna_QuadQuadGeo::Point(const Standard_Integer n) const 
2482 {
2483   if(!done)          {    StdFail_NotDone::Raise();        }
2484   if(n>nbint || n<1) {    Standard_DomainError::Raise();   }
2485   if(typeres==IntAna_PointAndCircle) {
2486     if(n!=1) { Standard_DomainError::Raise();  }
2487     if(param1==0.0) return(pt1);
2488     return(pt2);
2489   }
2490   else if(typeres==IntAna_Point) {
2491     if(n==1) return(pt1);
2492     return(pt2);
2493   }
2494
2495   // WNT (what can you expect from MicroSoft ?)
2496   return gp_Pnt(0,0,0);
2497 }
2498 //=======================================================================
2499 //function : Line
2500 //purpose  : Returns a Line
2501 //=======================================================================
2502   gp_Lin IntAna_QuadQuadGeo::Line(const Standard_Integer n) const 
2503 {
2504   if(!done)        {   StdFail_NotDone::Raise();   }
2505   if((n>nbint) || (n<1) || (typeres!=IntAna_Line)) {
2506     Standard_DomainError::Raise();
2507     }
2508   if(n==1) {  return(gp_Lin(pt1,dir1));   }
2509   else {      return(gp_Lin(pt2,dir2));   }
2510 }
2511 //=======================================================================
2512 //function : Circle
2513 //purpose  : Returns a Circle
2514 //=======================================================================
2515   gp_Circ IntAna_QuadQuadGeo::Circle(const Standard_Integer n) const 
2516 {
2517   if(!done) {    StdFail_NotDone::Raise();     }
2518   if(typeres==IntAna_PointAndCircle) {
2519     if(n!=1) { Standard_DomainError::Raise();  }
2520     if(param2==0.0) return(gp_Circ(DirToAx2(pt1,dir1),param1));
2521     return(gp_Circ(DirToAx2(pt2,dir2),param2));
2522   }
2523   else if((n>nbint) || (n<1) || (typeres!=IntAna_Circle)) {
2524     Standard_DomainError::Raise();
2525     }
2526   if      (n==1) { return(gp_Circ(DirToAx2(pt1,dir1),param1));}
2527   else if (n==2) { return(gp_Circ(DirToAx2(pt2,dir2),param2));}
2528   else if (n==3) { return(gp_Circ(DirToAx2(pt3,dir3),param3));}
2529   else           { return(gp_Circ(DirToAx2(pt4,dir4),param4));}
2530 }
2531
2532 //=======================================================================
2533 //function : Ellipse
2534 //purpose  : Returns a Elips  
2535 //=======================================================================
2536   gp_Elips IntAna_QuadQuadGeo::Ellipse(const Standard_Integer n) const
2537 {
2538   if(!done) {     StdFail_NotDone::Raise();     }
2539   if((n>nbint) || (n<1) || (typeres!=IntAna_Ellipse)) {
2540     Standard_DomainError::Raise();
2541   }
2542
2543   if(n==1) {
2544     Standard_Real R1=param1, R2=param1bis, aTmp;
2545     if (R1<R2) {
2546       aTmp=R1; R1=R2; R2=aTmp;
2547     }
2548     gp_Ax2 anAx2(pt1, dir1 ,dir2);
2549     gp_Elips anElips (anAx2, R1, R2);
2550     return anElips;
2551   }
2552   else {
2553     Standard_Real R1=param2, R2=param2bis, aTmp;
2554     if (R1<R2) {
2555       aTmp=R1; R1=R2; R2=aTmp;
2556     }
2557     gp_Ax2 anAx2(pt2, dir2 ,dir1);
2558     gp_Elips anElips (anAx2, R1, R2);
2559     return anElips;
2560   }
2561 }
2562 //=======================================================================
2563 //function : Parabola
2564 //purpose  : Returns a Parabola 
2565 //=======================================================================
2566   gp_Parab IntAna_QuadQuadGeo::Parabola(const Standard_Integer n) const 
2567 {
2568   if(!done) {
2569     StdFail_NotDone::Raise();
2570     }
2571   if (typeres!=IntAna_Parabola) {
2572     Standard_DomainError::Raise();
2573   }
2574   if((n>nbint) || (n!=1)) {
2575     Standard_OutOfRange::Raise();
2576   }
2577   return(gp_Parab(gp_Ax2( pt1
2578                          ,dir1
2579                          ,dir2)
2580                   ,param1));
2581 }
2582 //=======================================================================
2583 //function : Hyperbola
2584 //purpose  : Returns a Hyperbola  
2585 //=======================================================================
2586   gp_Hypr IntAna_QuadQuadGeo::Hyperbola(const Standard_Integer n) const 
2587 {
2588   if(!done) {
2589     StdFail_NotDone::Raise();
2590     }
2591   if((n>nbint) || (n<1) || (typeres!=IntAna_Hyperbola)) {
2592     Standard_DomainError::Raise();
2593     }
2594   if(n==1) {
2595     return(gp_Hypr(gp_Ax2( pt1
2596                           ,dir1
2597                           ,dir2)
2598                    ,param1,param1bis));
2599   }
2600   else {
2601     return(gp_Hypr(gp_Ax2( pt2
2602                           ,dir1
2603                           ,dir2.Reversed())
2604                    ,param2,param2bis));
2605   }
2606 }
2607 //=======================================================================
2608 //function : HasCommonGen
2609 //purpose  : 
2610 //=======================================================================
2611 Standard_Boolean IntAna_QuadQuadGeo::HasCommonGen() const
2612 {
2613   return myCommonGen;
2614 }
2615 //=======================================================================
2616 //function : PChar
2617 //purpose  : 
2618 //=======================================================================
2619 const gp_Pnt& IntAna_QuadQuadGeo::PChar() const
2620 {
2621   return myPChar;
2622 }
2623 //=======================================================================
2624 //function : RefineDir
2625 //purpose  : 
2626 //=======================================================================
2627 void RefineDir(gp_Dir& aDir)
2628 {
2629   Standard_Integer k, m, n;
2630   Standard_Real aC[3];
2631   //
2632   aDir.Coord(aC[0], aC[1], aC[2]);
2633   //
2634   m=0;
2635   n=0;
2636   for (k=0; k<3; ++k) {
2637     if (aC[k]==1. || aC[k]==-1.) {
2638       ++m;
2639     }
2640     else if (aC[k]!=0.) {
2641       ++n;
2642     }
2643   }
2644   //
2645   if (m && n) {
2646     Standard_Real aEps, aR1, aR2, aNum;
2647     //
2648     aEps=RealEpsilon();
2649     aR1=1.-aEps;
2650     aR2=1.+aEps;
2651     //
2652     for (k=0; k<3; ++k) {
2653       m=(aC[k]>0.);
2654       aNum=(m)? aC[k] : -aC[k];
2655       if (aNum>aR1 && aNum<aR2) {
2656         if (m) {
2657           aC[k]=1.;
2658         }          
2659         else {
2660           aC[k]=-1.;
2661         }
2662         //
2663         aC[(k+1)%3]=0.;
2664         aC[(k+2)%3]=0.;
2665         break;
2666       }
2667     }
2668     aDir.SetCoord(aC[0], aC[1], aC[2]);
2669   }
2670 }
2671
2672
2673