0024914: Micro edge is created during Boolean Operations
[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     if(DistA1A2<=Tol) {
960       if(RmR<=Tol) {
961         typeres=IntAna_Same;
962       }
963       else {
964         typeres=IntAna_Empty;
965       }
966     }
967     else {  //-- DistA1A2 > Tol
968       gp_Pnt P1=Cyl1.Location();
969       gp_Pnt P2t=Cyl2.Location();
970       gp_Pnt P2;
971       //-- P2t is projected on the plane (P1,DirCylX,DirCylY)
972       gp_Dir DirCyl = Cyl1.Position().Direction();
973       Standard_Real ProjP2OnDirCyl1=gp_Vec(DirCyl).Dot(gp_Vec(P1,P2t));
974       
975       P2.SetCoord( P2t.X() - ProjP2OnDirCyl1*DirCyl.X()
976                   ,P2t.Y() - ProjP2OnDirCyl1*DirCyl.Y()
977                   ,P2t.Z() - ProjP2OnDirCyl1*DirCyl.Z());
978       //-- 
979       Standard_Real R1pR2=R1+R2;
980       if(DistA1A2>(R1pR2+Tol)) { 
981         typeres=IntAna_Empty;
982         nbint=0;
983       }
984       else if(DistA1A2>(R1pR2)) {
985         //-- 1 Tangent line -------------------------------------OK
986         typeres=IntAna_Line;
987
988         nbint=1;
989         dir1=DirCyl;
990         Standard_Real R1_R1pR2=R1/R1pR2;
991         pt1.SetCoord( P1.X() + R1_R1pR2 * (P2.X()-P1.X())
992                      ,P1.Y() + R1_R1pR2 * (P2.Y()-P1.Y())
993                      ,P1.Z() + R1_R1pR2 * (P2.Z()-P1.Z()));
994         
995       }
996       else if(DistA1A2>RmR) {
997         //-- 2 lines ---------------------------------------------OK
998         typeres=IntAna_Line;
999         nbint=2;
1000         dir1=DirCyl;
1001         gp_Vec P1P2(P1,P2);
1002         gp_Dir DirA1A2=gp_Dir(P1P2);
1003         gp_Dir Ortho_dir1_P1P2 = dir1.Crossed(DirA1A2);
1004         dir2=dir1;
1005         Standard_Real Alpha=0.5*(R1*R1-R2*R2+DistA1A2*DistA1A2)/(DistA1A2);       
1006
1007 //         Standard_Real Beta = Sqrt(R1*R1-Alpha*Alpha);
1008         Standard_Real anSqrtArg = R1*R1-Alpha*Alpha;
1009         Standard_Real Beta = (anSqrtArg > 0.) ? Sqrt(anSqrtArg) : 0.;
1010         
1011         if((Beta+Beta)<Tol) { 
1012           nbint=1;
1013           pt1.SetCoord( P1.X() + Alpha*DirA1A2.X()
1014                        ,P1.Y() + Alpha*DirA1A2.Y()
1015                        ,P1.Z() + Alpha*DirA1A2.Z());
1016         }
1017         else { 
1018           pt1.SetCoord( P1.X() + Alpha*DirA1A2.X() + Beta*Ortho_dir1_P1P2.X()
1019                        ,P1.Y() + Alpha*DirA1A2.Y() + Beta*Ortho_dir1_P1P2.Y()
1020                        ,P1.Z() + Alpha*DirA1A2.Z() + Beta*Ortho_dir1_P1P2.Z() );
1021           
1022           pt2.SetCoord( P1.X() + Alpha*DirA1A2.X() - Beta*Ortho_dir1_P1P2.X()
1023                        ,P1.Y() + Alpha*DirA1A2.Y() - Beta*Ortho_dir1_P1P2.Y()
1024                        ,P1.Z() + Alpha*DirA1A2.Z() - Beta*Ortho_dir1_P1P2.Z());
1025         }
1026       }
1027       else if(DistA1A2>(RmR-Tol)) {
1028         //-- 1 Tangent ------------------------------------------OK
1029         typeres=IntAna_Line;
1030         nbint=1;
1031         dir1=DirCyl;
1032         Standard_Real R1_RmR=R1/RmR;
1033
1034         if(R1 < R2) R1_RmR = -R1_RmR;
1035
1036         pt1.SetCoord( P1.X() + R1_RmR * (P2.X()-P1.X())
1037                      ,P1.Y() + R1_RmR * (P2.Y()-P1.Y())
1038                      ,P1.Z() + R1_RmR * (P2.Z()-P1.Z()));
1039       }
1040       else {
1041         nbint=0;
1042         typeres=IntAna_Empty;
1043       }
1044     }
1045   }
1046   else { //-- No Parallel Axis ---------------------------------OK
1047     if((RmR_Relative<=myEPSILON_CYLINDER_DELTA_RADIUS) 
1048        && (DistA1A2 <= myEPSILON_CYLINDER_DELTA_DISTANCE)) {
1049       //-- PI/2 between the two axis   and   Intersection  
1050       //-- and identical radius
1051       typeres=IntAna_Ellipse;
1052       nbint=2;
1053       gp_Dir DirCyl1=Cyl1.Position().Direction();
1054       gp_Dir DirCyl2=Cyl2.Position().Direction();
1055       pt1=pt2=A1A2.PtIntersect();
1056       
1057       Standard_Real A=DirCyl1.Angle(DirCyl2);
1058       Standard_Real B;
1059       B=Abs(Sin(0.5*(M_PI-A)));
1060       A=Abs(Sin(0.5*A));
1061       
1062       if(A==0.0 || B==0.0) {
1063         typeres=IntAna_Same;
1064         return;
1065       }
1066       
1067       
1068       gp_Vec dircyl1(DirCyl1);gp_Vec dircyl2(DirCyl2);
1069       dir1 = gp_Dir(dircyl1.Added(dircyl2));
1070       dir2 = gp_Dir(dircyl1.Subtracted(dircyl2));
1071         
1072       param2   = Cyl1.Radius() / A;
1073       param1   = Cyl1.Radius() / B;
1074       param2bis= param1bis = Cyl1.Radius();
1075       if(param1 < param1bis) {
1076         A=param1; param1=param1bis; param1bis=A;
1077       }
1078       if(param2 < param2bis) {
1079         A=param2; param2=param2bis; param2bis=A;
1080       }
1081     }
1082     else {
1083       if(Abs(DistA1A2-Cyl1.Radius()-Cyl2.Radius())<Tol) { 
1084         typeres = IntAna_Point;
1085         Standard_Real d,p1,p2;
1086
1087         gp_Dir D1 = Cyl1.Axis().Direction();
1088         gp_Dir D2 = Cyl2.Axis().Direction();
1089         A1A2.Distance(d,p1,p2);
1090         gp_Pnt P = Cyl1.Axis().Location();
1091         gp_Pnt P1(P.X() - p1*D1.X(),
1092                   P.Y() - p1*D1.Y(),
1093                   P.Z() - p1*D1.Z());
1094         P = Cyl2.Axis().Location();
1095         gp_Pnt P2(P.X() - p2*D2.X(),
1096                   P.Y() - p2*D2.Y(),
1097                   P.Z() - p2*D2.Z());
1098         gp_Vec P1P2(P1,P2);
1099         D1=gp_Dir(P1P2);
1100         p1=Cyl1.Radius();
1101         pt1.SetCoord(P1.X() + p1*D1.X(),
1102                      P1.Y() + p1*D1.Y(),
1103                      P1.Z() + p1*D1.Z());
1104         nbint = 1;
1105       }
1106       else {
1107         typeres=IntAna_NoGeometricSolution;
1108       }
1109     }
1110   }
1111 }
1112 //=======================================================================
1113 //function : IntAna_QuadQuadGeo
1114 //purpose  : Cylinder - Cone
1115 //=======================================================================
1116 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
1117                                        const gp_Cone& Con,
1118                                        const Standard_Real Tol) 
1119 : done(Standard_False),
1120   nbint(0),
1121   typeres(IntAna_Empty),
1122   pt1(0,0,0),
1123   pt2(0,0,0),
1124   pt3(0,0,0),
1125   pt4(0,0,0),
1126   param1(0),
1127   param2(0),
1128   param3(0),
1129   param4(0),
1130   param1bis(0),
1131   param2bis(0),
1132   myCommonGen(Standard_False),
1133   myPChar(0,0,0)
1134 {
1135   InitTolerances();
1136   Perform(Cyl,Con,Tol);
1137 }
1138 //=======================================================================
1139 //function : Perform
1140 //purpose  : 
1141 //=======================================================================
1142   void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl,
1143                                    const gp_Cone& Con,
1144                                    const Standard_Real ) 
1145 {
1146   done=Standard_True;
1147   AxeOperator A1A2(Cyl.Axis(),Con.Axis());
1148   if(A1A2.Same()) {
1149     gp_Pnt Pt=Con.Apex();
1150     Standard_Real dist=Cyl.Radius()/(Tan(Con.SemiAngle()));
1151     gp_Dir dir=Cyl.Position().Direction();
1152     pt1.SetCoord( Pt.X() + dist*dir.X()
1153                  ,Pt.Y() + dist*dir.Y()
1154                  ,Pt.Z() + dist*dir.Z());
1155     pt2.SetCoord( Pt.X() - dist*dir.X()
1156                  ,Pt.Y() - dist*dir.Y()
1157                  ,Pt.Z() - dist*dir.Z());
1158     dir1=dir2=dir;
1159     param1=param2=Cyl.Radius();
1160     nbint=2;
1161     typeres=IntAna_Circle;
1162
1163   }
1164   else {
1165     typeres=IntAna_NoGeometricSolution;
1166   }
1167 }
1168 //=======================================================================
1169 //function : 
1170 //purpose  : Cylinder - Sphere
1171 //=======================================================================
1172   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
1173                                          const gp_Sphere& Sph,
1174                                          const Standard_Real Tol) 
1175 : done(Standard_False),
1176   nbint(0),
1177   typeres(IntAna_Empty),
1178   pt1(0,0,0),
1179   pt2(0,0,0),
1180   pt3(0,0,0),
1181   pt4(0,0,0),
1182   param1(0),
1183   param2(0),
1184   param3(0),
1185   param4(0),
1186   param1bis(0),
1187   param2bis(0),
1188   myCommonGen(Standard_False),
1189   myPChar(0,0,0)
1190 {
1191   InitTolerances();
1192   Perform(Cyl,Sph,Tol);
1193 }
1194 //=======================================================================
1195 //function : Perform
1196 //purpose  : 
1197 //=======================================================================
1198   void IntAna_QuadQuadGeo::Perform( const gp_Cylinder& Cyl
1199                                    ,const gp_Sphere& Sph
1200                                    ,const Standard_Real)
1201 {
1202   done=Standard_True;
1203   gp_Pnt Pt=Sph.Location();
1204   AxeOperator A1A2(Cyl.Axis(),Sph.Position().Axis());
1205   if((A1A2.Intersect()  && Pt.Distance(A1A2.PtIntersect())==0.0 )
1206      || (A1A2.Same()))      {
1207     if(Sph.Radius() < Cyl.Radius()) { 
1208       typeres = IntAna_Empty;
1209     }
1210     else { 
1211       Standard_Real dist=Sqrt( Sph.Radius() * Sph.Radius() - Cyl.Radius() * Cyl.Radius() );
1212       gp_Dir dir=Cyl.Position().Direction();
1213       dir1 = dir2 = dir;
1214       typeres=IntAna_Circle;
1215       pt1.SetCoord( Pt.X() + dist*dir.X()
1216                    ,Pt.Y() + dist*dir.Y()
1217                    ,Pt.Z() + dist*dir.Z());
1218       nbint=1;
1219       param1 = Cyl.Radius();
1220       if(dist>RealEpsilon()) {
1221         pt2.SetCoord( Pt.X() - dist*dir.X()
1222                      ,Pt.Y() - dist*dir.Y()
1223                      ,Pt.Z() - dist*dir.Z());
1224         param2=Cyl.Radius();
1225         nbint=2;
1226       }
1227     }
1228   }
1229   else {
1230     typeres=IntAna_NoGeometricSolution; 
1231   }
1232 }
1233
1234 //=======================================================================
1235 //function : IntAna_QuadQuadGeo
1236 //purpose  : Cone - Cone
1237 //=======================================================================
1238   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cone& Con1,
1239                                          const gp_Cone& Con2,
1240                                          const Standard_Real Tol) 
1241 : done(Standard_False),
1242   nbint(0),
1243   typeres(IntAna_Empty),
1244   pt1(0,0,0),
1245   pt2(0,0,0),
1246   pt3(0,0,0),
1247   pt4(0,0,0),
1248   param1(0),
1249   param2(0),
1250   param3(0),
1251   param4(0),
1252   param1bis(0),
1253   param2bis(0),
1254   myCommonGen(Standard_False),
1255   myPChar(0,0,0)
1256 {
1257   InitTolerances();
1258   Perform(Con1,Con2,Tol);
1259 }
1260 //
1261 //=======================================================================
1262 //function : Perform
1263 //purpose  : 
1264 //=======================================================================
1265   void IntAna_QuadQuadGeo::Perform(const gp_Cone& Con1,
1266                                    const gp_Cone& Con2,
1267                                    const Standard_Real Tol) 
1268 {
1269   done=Standard_True;
1270   //
1271   Standard_Real tg1, tg2, aDA1A2, aTol2;
1272   gp_Pnt aPApex1, aPApex2;
1273
1274   Standard_Real TOL_APEX_CONF = 1.e-10;
1275   
1276   //
1277   tg1=Tan(Con1.SemiAngle());
1278   tg2=Tan(Con2.SemiAngle());
1279
1280   if((tg1 * tg2) < 0.) {
1281     tg2 = -tg2;
1282   }
1283   //
1284   aTol2=Tol*Tol;
1285   aPApex1=Con1.Apex();
1286   aPApex2=Con2.Apex();
1287   aDA1A2=aPApex1.SquareDistance(aPApex2);
1288   //
1289   AxeOperator A1A2(Con1.Axis(),Con2.Axis());
1290   //
1291   // 1
1292   if(A1A2.Same()) {
1293     //-- two circles 
1294     Standard_Real x;
1295     gp_Pnt P=Con1.Apex();
1296     gp_Dir D=Con1.Position().Direction();
1297     Standard_Real d=gp_Vec(D).Dot(gp_Vec(P,Con2.Apex()));
1298     
1299     if(Abs(tg1-tg2)>myEPSILON_ANGLE_CONE) { 
1300       if (fabs(d) < TOL_APEX_CONF) {
1301         typeres = IntAna_Point;
1302         nbint = 1;
1303         pt1 = P;
1304         return;
1305       }
1306       x=(d*tg2)/(tg1+tg2);
1307       pt1.SetCoord( P.X() + x*D.X()
1308                    ,P.Y() + x*D.Y()
1309                    ,P.Z() + x*D.Z());
1310       param1=Abs(x*tg1);
1311
1312       x=(d*tg2)/(tg2-tg1);
1313       pt2.SetCoord( P.X() + x*D.X()
1314                    ,P.Y() + x*D.Y()
1315                    ,P.Z() + x*D.Z());
1316       param2=Abs(x*tg1);
1317       dir1 = dir2 = D;
1318       nbint=2;
1319       typeres=IntAna_Circle;
1320     }
1321     else {
1322       if (fabs(d) < TOL_APEX_CONF) { 
1323         typeres=IntAna_Same;
1324       }
1325       else {
1326         typeres=IntAna_Circle;
1327         nbint=1;
1328         x=d*0.5;
1329         pt1.SetCoord( P.X() + x*D.X()
1330                      ,P.Y() + x*D.Y()
1331                      ,P.Z() + x*D.Z());
1332         param1 = Abs(x * tg1);
1333         dir1 = D;
1334       }
1335     }
1336   } //-- fin A1A2.Same
1337   // 2
1338   else if((Abs(tg1-tg2)<myEPSILON_ANGLE_CONE) && (A1A2.Parallel())) {
1339     //-- voir AnVer12mai98
1340     Standard_Real DistA1A2=A1A2.Distance();
1341     gp_Dir DA1=Con1.Position().Direction();
1342     gp_Vec O1O2(Con1.Apex(),Con2.Apex());
1343     gp_Dir O1O2n(O1O2); // normalization of the vector before projection
1344     Standard_Real O1O2_DA1=gp_Vec(DA1).Dot(gp_Vec(O1O2n));
1345
1346     gp_Vec O1_Proj_A2(O1O2n.X()-O1O2_DA1*DA1.X(),
1347                       O1O2n.Y()-O1O2_DA1*DA1.Y(),
1348                       O1O2n.Z()-O1O2_DA1*DA1.Z());
1349     gp_Dir DB1=gp_Dir(O1_Proj_A2);
1350
1351     Standard_Real yO1O2=O1O2.Dot(gp_Vec(DA1));
1352     Standard_Real ABSTG1 = Abs(tg1);
1353     Standard_Real X2 = (DistA1A2/ABSTG1 - yO1O2)*0.5;
1354     Standard_Real X1 = X2+yO1O2;
1355     
1356     gp_Pnt P1(Con1.Apex().X() + X1*( DA1.X() + ABSTG1*DB1.X()),
1357               Con1.Apex().Y() + X1*( DA1.Y() + ABSTG1*DB1.Y()), 
1358               Con1.Apex().Z() + X1*( DA1.Z() + ABSTG1*DB1.Z()));
1359
1360     gp_Pnt MO1O2(0.5*(Con1.Apex().X()+Con2.Apex().X()),
1361                  0.5*(Con1.Apex().Y()+Con2.Apex().Y()),
1362                  0.5*(Con1.Apex().Z()+Con2.Apex().Z()));
1363     gp_Vec P1MO1O2(P1,MO1O2);
1364     
1365     gp_Dir DA1_X_DB1=DA1.Crossed(DB1);
1366     gp_Dir OrthoPln =  DA1_X_DB1.Crossed(gp_Dir(P1MO1O2));
1367     
1368     IntAna_QuadQuadGeo INTER_QUAD_PLN(gp_Pln(P1,OrthoPln),Con1,Tol,Tol);
1369       if(INTER_QUAD_PLN.IsDone()) {
1370       switch(INTER_QUAD_PLN.TypeInter()) {
1371       case IntAna_Ellipse:         {
1372         typeres=IntAna_Ellipse;
1373         gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
1374         pt1 = E.Location();
1375         dir1 = E.Position().Direction();
1376         dir2 = E.Position().XDirection();
1377         param1 = E.MajorRadius();
1378         param1bis = E.MinorRadius();
1379         nbint = 1;
1380         break; 
1381       }
1382       case IntAna_Circle: {
1383         typeres=IntAna_Circle;
1384         gp_Circ C=INTER_QUAD_PLN.Circle(1);
1385         pt1 = C.Location();
1386         dir1 = C.Position().XDirection();
1387         dir2 = C.Position().YDirection();
1388         param1 = C.Radius();
1389         nbint = 1;
1390         break;
1391       }
1392       case IntAna_Hyperbola: {
1393         typeres=IntAna_Hyperbola;
1394         gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
1395         pt1 = pt2 = H.Location();
1396         dir1 = H.Position().Direction();
1397         dir2 = H.Position().XDirection();
1398         param1 = param2 = H.MajorRadius();
1399         param1bis = param2bis = H.MinorRadius();
1400         nbint = 2;
1401         break;
1402       }
1403       case IntAna_Line: {
1404         typeres=IntAna_Line;
1405         gp_Lin H=INTER_QUAD_PLN.Line(1);
1406         pt1 = pt2 = H.Location();
1407         dir1 = dir2 = H.Position().Direction();
1408         param1 = param2 = 0.0;
1409         param1bis = param2bis = 0.0;
1410         nbint = 2;
1411         break;
1412       }
1413       default:
1414         typeres=IntAna_NoGeometricSolution; 
1415       }
1416     }
1417   }// else if((Abs(tg1-tg2)<EPSILON_ANGLE_CONE) && (A1A2.Parallel()))
1418   // 3
1419   else if (aDA1A2<aTol2) {
1420     //
1421     // When apices are coinsided there can be 3 possible cases
1422     // 3.1 - empty solution (iRet=0)
1423     // 3.2 - one line  when cone1 touches cone2 (iRet=1)
1424     // 3.3 - two lines when cone1 intersects cone2 (iRet=2)
1425     //
1426     Standard_Integer iRet;
1427     Standard_Real aGamma, aBeta1, aBeta2;
1428     Standard_Real aD1, aR1, aTgBeta1, aTgBeta2, aHalfPI;
1429     Standard_Real aCosGamma, aSinGamma, aDx, aR2, aRD2, aD2;
1430     gp_Pnt2d aP0, aPA1, aP1, aPA2;
1431     gp_Vec2d aVAx2;
1432     gp_Ax1 aAx1, aAx2;
1433     //
1434     // Preliminary analysis. Determination of iRet
1435     //
1436     iRet=0;
1437     aHalfPI=0.5*M_PI;
1438     aD1=1.;
1439     aPA1.SetCoord(aD1, 0.);
1440     aP0.SetCoord(0., 0.);
1441     //
1442     aAx1=Con1.Axis();
1443     aAx2=Con2.Axis();
1444     aGamma=aAx1.Angle(aAx2);
1445     if (aGamma>aHalfPI){
1446       aGamma=M_PI-aGamma;
1447     }
1448     aCosGamma=Cos(aGamma);
1449     aSinGamma=Sin(aGamma);
1450     //
1451     aBeta1=Con1.SemiAngle();
1452     aTgBeta1=Tan(aBeta1);
1453     aTgBeta1=Abs(aTgBeta1);
1454     //
1455     aBeta2=Con2.SemiAngle();
1456     aTgBeta2=Tan(aBeta2);
1457     aTgBeta2=Abs(aTgBeta2);
1458     //
1459     aR1=aD1*aTgBeta1;
1460     aP1.SetCoord(aD1, aR1);
1461     //
1462     // PA2
1463     aVAx2.SetCoord(aCosGamma, aSinGamma);
1464     gp_Dir2d aDAx2(aVAx2);
1465     gp_Lin2d aLAx2(aP0, aDAx2);
1466     //
1467     gp_Vec2d aV(aP0, aP1);
1468     aDx=aVAx2.Dot(aV);
1469     aPA2=aP0.Translated(aDx*aDAx2);
1470     //
1471     // aR2
1472     aDx=aPA2.Distance(aP0);
1473     aR2=aDx*aTgBeta2;
1474     //
1475     // aRD2
1476     aRD2=aPA2.Distance(aP1);
1477     //
1478     if (aRD2>(aR2+Tol)) {
1479       iRet=0;
1480       typeres=IntAna_Empty; //nothing
1481       return;
1482     }
1483     //
1484     iRet=1; //touch case => 1 line
1485     if (aRD2<(aR2-Tol)) {
1486       iRet=2;//intersection => couple of lines
1487     }
1488     //
1489     // Finding the solution in 3D
1490     //
1491     Standard_Real aDa;
1492     gp_Pnt aQApex1, aQA1, aQA2, aQX, aQX1, aQX2;
1493     gp_Dir aD3Ax1, aD3Ax2;
1494     gp_Lin aLin;
1495     IntAna_QuadQuadGeo aIntr;
1496     //
1497     aQApex1=Con1.Apex();
1498     aD3Ax1=aAx1.Direction(); 
1499     aQA1.SetCoord(aQApex1.X()+aD1*aD3Ax1.X(),
1500                   aQApex1.Y()+aD1*aD3Ax1.Y(),
1501                   aQApex1.Z()+aD1*aD3Ax1.Z());
1502     //
1503     aDx=aD3Ax1.Dot(aAx2.Direction());
1504     if (aDx<0.) {
1505       aAx2.Reverse();
1506     }
1507     aD3Ax2=aAx2.Direction();
1508     //
1509     aD2=aD1*sqrt((1.+aTgBeta1*aTgBeta1)/(1.+aTgBeta2*aTgBeta2));
1510     //
1511     aQA2.SetCoord(aQApex1.X()+aD2*aD3Ax2.X(),
1512                   aQApex1.Y()+aD2*aD3Ax2.Y(),
1513                   aQApex1.Z()+aD2*aD3Ax2.Z());
1514     //
1515     gp_Pln aPln1(aQA1, aD3Ax1);
1516     gp_Pln aPln2(aQA2, aD3Ax2);
1517     //
1518     aIntr.Perform(aPln1, aPln2, Tol, Tol);
1519     if (!aIntr.IsDone()) {
1520       iRet=-1; // just in case. it must not be so
1521       typeres=IntAna_NoGeometricSolution; 
1522       return;
1523     }
1524     //
1525     aLin=aIntr.Line(1);
1526     const gp_Dir& aDLin=aLin.Direction();
1527     gp_Vec aVLin(aDLin);
1528     gp_Pnt aOrig=aLin.Location();
1529     gp_Vec aVr(aQA1, aOrig);
1530     aDx=aVLin.Dot(aVr);
1531     aQX=aOrig.Translated(aDx*aVLin);
1532     //
1533     // Final part
1534     //
1535     typeres=IntAna_Line; 
1536     //
1537     param1=0.;
1538     param2 =0.;
1539     param1bis=0.;
1540     param2bis=0.;
1541     //
1542     if (iRet==1) {
1543       // one line
1544       nbint=1;
1545       pt1=aQApex1;
1546       gp_Vec aVX(aQApex1, aQX);
1547       dir1=gp_Dir(aVX);
1548     }
1549     
1550     else {//iRet=2 
1551       // two lines
1552       nbint=2;
1553       aDa=aQA1.Distance(aQX);
1554       aDx=sqrt(aR1*aR1-aDa*aDa);
1555       aQX1=aQX.Translated(aDx*aVLin);
1556       aQX2=aQX.Translated(-aDx*aVLin);
1557       //
1558       pt1=aQApex1;
1559       pt2=aQApex1;
1560       gp_Vec aVX1(aQApex1, aQX1);
1561       dir1=gp_Dir(aVX1);
1562       gp_Vec aVX2(aQApex1, aQX2);
1563       dir2=gp_Dir(aVX2);
1564     }
1565   } //else if (aDA1A2<aTol2) {
1566   //Case when cones have common generatrix
1567   else if(A1A2.Intersect()) {
1568     //Check if apex of one cone belongs another one
1569     Standard_Real u, v, tol2 = Tol*Tol;
1570     ElSLib::Parameters(Con2, aPApex1, u, v);
1571     gp_Pnt p = ElSLib::Value(u, v, Con2);
1572     if(aPApex1.SquareDistance(p) > tol2) {
1573       typeres=IntAna_NoGeometricSolution; 
1574       return;
1575     }
1576     //
1577     ElSLib::Parameters(Con1, aPApex2, u, v);
1578     p = ElSLib::Value(u, v, Con1);
1579     if(aPApex2.SquareDistance(p) > tol2) {
1580       typeres=IntAna_NoGeometricSolution; 
1581       return;
1582     }
1583
1584     //Cones have a common generatrix passing through apexes
1585     myCommonGen = Standard_True;
1586
1587     //common generatrix of cones
1588     gp_Lin aGen(aPApex1, gp_Dir(gp_Vec(aPApex1, aPApex2)));
1589
1590     //Intersection point of axes
1591     gp_Pnt aPAxeInt = A1A2.PtIntersect();
1592
1593     //Characteristic point of intersection curve
1594     u = ElCLib::Parameter(aGen, aPAxeInt);
1595     myPChar = ElCLib::Value(u, aGen);
1596
1597
1598     //Other generatrixes of cones laying in maximal plane
1599     gp_Lin aGen1 = aGen.Rotated(Con1.Axis(), M_PI); 
1600     gp_Lin aGen2 = aGen.Rotated(Con2.Axis(), M_PI); 
1601     //
1602     //Intersection point of generatrixes
1603     gp_Dir aN; //solution plane normal
1604     gp_Dir aD1 = aGen1.Direction();
1605
1606     gp_Dir aD2(aD1.Crossed(aGen.Direction()));
1607
1608     if(aD1.IsParallel(aGen2.Direction(), Precision::Angular())) {
1609       aN = aD1.Crossed(aD2);
1610     }
1611     else if(aGen1.SquareDistance(aGen2) > tol2) {
1612       //Something wrong ???
1613       typeres=IntAna_NoGeometricSolution; 
1614       return;
1615     }
1616     else {
1617       gp_Dir D1 = aGen1.Position().Direction();
1618       gp_Dir D2 = aGen2.Position().Direction();
1619       gp_Pnt O1 = aGen1.Location();
1620       gp_Pnt O2 = aGen2.Location();
1621       Standard_Real D1DotD2 = D1.Dot(D2);
1622       Standard_Real aSin = 1.-D1DotD2*D1DotD2;
1623       gp_Vec O1O2 (O1,O2);
1624       Standard_Real U2 = (D1.XYZ()*(O1O2.Dot(D1))-(O1O2.XYZ())).Dot(D2.XYZ());
1625       U2 /= aSin;
1626       gp_Pnt aPGint(ElCLib::Value(U2, aGen2));
1627     
1628       aD1 = gp_Dir(gp_Vec(aPGint, myPChar));
1629       aN = aD1.Crossed(aD2);
1630     }
1631     //Plane that must contain intersection curves
1632     gp_Pln anIntPln(myPChar, aN);
1633
1634     IntAna_QuadQuadGeo INTER_QUAD_PLN(anIntPln,Con1,Tol,Tol);
1635
1636       if(INTER_QUAD_PLN.IsDone()) {
1637       switch(INTER_QUAD_PLN.TypeInter()) {
1638       case IntAna_Ellipse:         {
1639         typeres=IntAna_Ellipse;
1640         gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
1641         pt1 = E.Location();
1642         dir1 = E.Position().Direction();
1643         dir2 = E.Position().XDirection();
1644         param1 = E.MajorRadius();
1645         param1bis = E.MinorRadius();
1646         nbint = 1;
1647         break; 
1648       }
1649       case IntAna_Circle: {
1650         typeres=IntAna_Circle;
1651         gp_Circ C=INTER_QUAD_PLN.Circle(1);
1652         pt1 = C.Location();
1653         dir1 = C.Position().XDirection();
1654         dir2 = C.Position().YDirection();
1655         param1 = C.Radius();
1656         nbint = 1;
1657         break;
1658       }
1659       case IntAna_Parabola: {
1660         typeres=IntAna_Parabola;
1661         gp_Parab Prb=INTER_QUAD_PLN.Parabola(1);
1662         pt1 = Prb.Location();
1663         dir1 = Prb.Position().Direction();
1664         dir2 = Prb.Position().XDirection();
1665         param1 = Prb.Focal();
1666         nbint = 1;
1667         break;
1668       }
1669       case IntAna_Hyperbola: {
1670         typeres=IntAna_Hyperbola;
1671         gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
1672         pt1 = pt2 = H.Location();
1673         dir1 = H.Position().Direction();
1674         dir2 = H.Position().XDirection();
1675         param1 = param2 = H.MajorRadius();
1676         param1bis = param2bis = H.MinorRadius();
1677         nbint = 2;
1678         break;
1679       }
1680       default:
1681         typeres=IntAna_NoGeometricSolution; 
1682       }
1683     }
1684   }
1685   
1686   else {
1687     typeres=IntAna_NoGeometricSolution; 
1688   }
1689 }
1690 //=======================================================================
1691 //function : IntAna_QuadQuadGeo
1692 //purpose  : Sphere - Cone
1693 //=======================================================================
1694   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Sphere& Sph,
1695                                          const gp_Cone& Con,
1696                                          const Standard_Real Tol) 
1697 : done(Standard_False),
1698   nbint(0),
1699   typeres(IntAna_Empty),
1700   pt1(0,0,0),
1701   pt2(0,0,0),
1702   pt3(0,0,0),
1703   pt4(0,0,0),
1704   param1(0),
1705   param2(0),
1706   param3(0),
1707   param4(0),
1708   param1bis(0),
1709   param2bis(0),
1710   myCommonGen(Standard_False),
1711   myPChar(0,0,0)
1712 {
1713   InitTolerances();
1714   Perform(Sph,Con,Tol);
1715 }
1716 //=======================================================================
1717 //function : Perform
1718 //purpose  : 
1719 //=======================================================================
1720   void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph,
1721                                    const gp_Cone& Con,
1722                                    const Standard_Real)
1723 {
1724   
1725   //
1726   done=Standard_True;
1727   //
1728   AxeOperator A1A2(Con.Axis(),Sph.Position().Axis());
1729   gp_Pnt Pt=Sph.Location();
1730   //
1731   if((A1A2.Intersect() && (Pt.Distance(A1A2.PtIntersect())==0.0))
1732      || A1A2.Same()) {
1733     gp_Pnt ConApex= Con.Apex();
1734     Standard_Real dApexSphCenter=Pt.Distance(ConApex); 
1735     gp_Dir ConDir;
1736     if(dApexSphCenter>RealEpsilon()) { 
1737       ConDir = gp_Dir(gp_Vec(ConApex,Pt));
1738     }
1739     else { 
1740       ConDir = Con.Position().Direction();
1741     }
1742     
1743     Standard_Real Rad=Sph.Radius();
1744     Standard_Real tga=Tan(Con.SemiAngle());
1745
1746
1747     //-- 2 circles
1748     //-- x: Roots of    (x**2 + y**2 = Rad**2)
1749     //--                tga = y / (x+dApexSphCenter)
1750     Standard_Real tgatga = tga * tga;
1751     math_DirectPolynomialRoots Eq( 1.0+tgatga
1752                                   ,2.0*tgatga*dApexSphCenter
1753                                   ,-Rad*Rad + dApexSphCenter*dApexSphCenter*tgatga);
1754     if(Eq.IsDone()) {
1755       Standard_Integer nbsol=Eq.NbSolutions();
1756       if(nbsol==0) {
1757         typeres=IntAna_Empty;
1758       }
1759       else { 
1760         typeres=IntAna_Circle;
1761         if(nbsol>=1) {
1762           Standard_Real x = Eq.Value(1);
1763           Standard_Real dApexSphCenterpx = dApexSphCenter+x;
1764           nbint=1;
1765           pt1.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
1766                        ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
1767                        ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
1768           param1 = tga * dApexSphCenterpx;
1769           param1 = Abs(param1);
1770           dir1 = ConDir;
1771           if(param1<=myEPSILON_MINI_CIRCLE_RADIUS) {
1772             typeres=IntAna_PointAndCircle;
1773             param1=0.0;
1774           }
1775         }
1776         if(nbsol>=2) {
1777           Standard_Real x=Eq.Value(2);
1778           Standard_Real dApexSphCenterpx = dApexSphCenter+x;
1779           nbint=2;
1780           pt2.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
1781                        ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
1782                        ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
1783           param2 = tga * dApexSphCenterpx;
1784           param2 = Abs(param2);
1785           dir2=ConDir;
1786           if(param2<=myEPSILON_MINI_CIRCLE_RADIUS) {
1787             typeres=IntAna_PointAndCircle;
1788             param2=0.0;
1789           }
1790         }
1791       }
1792     }
1793     else {
1794       done=Standard_False;
1795     }
1796   }
1797   else {
1798     typeres=IntAna_NoGeometricSolution; 
1799   }
1800 }
1801
1802 //=======================================================================
1803 //function : IntAna_QuadQuadGeo
1804 //purpose  : Sphere - Sphere
1805 //=======================================================================
1806   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo( const gp_Sphere& Sph1
1807                                          ,const gp_Sphere& Sph2
1808                                          ,const Standard_Real Tol) 
1809 : done(Standard_False),
1810   nbint(0),
1811   typeres(IntAna_Empty),
1812   pt1(0,0,0),
1813   pt2(0,0,0),
1814   pt3(0,0,0),
1815   pt4(0,0,0),
1816   param1(0),
1817   param2(0),
1818   param3(0),
1819   param4(0),
1820   param1bis(0),
1821   param2bis(0),
1822   myCommonGen(Standard_False),
1823   myPChar(0,0,0)
1824 {
1825   InitTolerances();
1826   Perform(Sph1,Sph2,Tol);
1827 }
1828 //=======================================================================
1829 //function : Perform
1830 //purpose  : 
1831 //=======================================================================
1832   void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph1,
1833                                    const gp_Sphere& Sph2,
1834                                    const Standard_Real Tol)   
1835 {
1836   done=Standard_True;
1837   gp_Pnt O1=Sph1.Location();
1838   gp_Pnt O2=Sph2.Location();
1839   Standard_Real dO1O2=O1.Distance(O2);
1840   Standard_Real R1=Sph1.Radius();
1841   Standard_Real R2=Sph2.Radius();
1842   Standard_Real Rmin,Rmax;
1843   typeres=IntAna_Empty;
1844   param2bis=0.0; //-- pour eviter param2bis not used .... 
1845
1846   if(R1>R2) { Rmin=R2; Rmax=R1; } else { Rmin=R1; Rmax=R2; }
1847   
1848   if(dO1O2<=Tol && (Abs(R1-R2) <= Tol)) {
1849     typeres = IntAna_Same;
1850   }
1851   else { 
1852     if(dO1O2<=Tol) { return; } 
1853     gp_Dir Dir=gp_Dir(gp_Vec(O1,O2));
1854     Standard_Real t = Rmax - dO1O2 - Rmin;
1855
1856     //----------------------------------------------------------------------
1857     //--        |----------------- R1 --------------------|
1858     //--        |----dO1O2-----|-----------R2----------|
1859     //--                                            --->--<-- t
1860     //--
1861     //--        |------ R1 ------|---------dO1O2----------|
1862     //--     |-------------------R2-----------------------|
1863     //--  --->--<-- t
1864     //----------------------------------------------------------------------
1865     if(t >= 0.0  && t <=Tol) { 
1866       typeres = IntAna_Point;
1867       nbint = 1;
1868       Standard_Real t2;
1869       if(R1==Rmax) t2=(R1 + (R2 + dO1O2)) * 0.5;
1870       else         t2=(-R1+(dO1O2-R2))*0.5;
1871         
1872       pt1.SetCoord( O1.X() + t2*Dir.X()
1873                    ,O1.Y() + t2*Dir.Y()
1874                    ,O1.Z() + t2*Dir.Z());
1875     }
1876     else  {
1877       //-----------------------------------------------------------------
1878       //--        |----------------- dO1O2 --------------------|
1879       //--        |----R1-----|-----------R2----------|-Tol-|
1880       //--                                            
1881       //--        |----------------- Rmax --------------------|
1882       //--        |----Rmin----|-------dO1O2-------|-Tol-|
1883       //--                                            
1884       //-----------------------------------------------------------------
1885       if((dO1O2 > (R1+R2+Tol)) || (Rmax > (dO1O2+Rmin+Tol))) {
1886         typeres=IntAna_Empty;
1887       }
1888       else {
1889         //---------------------------------------------------------------
1890         //--     
1891         //--
1892         //---------------------------------------------------------------
1893         Standard_Real Alpha=0.5*(R1*R1-R2*R2+dO1O2*dO1O2)/(dO1O2);       
1894         Standard_Real Beta = R1*R1-Alpha*Alpha;
1895         Beta = (Beta>0.0)? Sqrt(Beta) : 0.0;
1896         
1897         if(Beta<= myEPSILON_MINI_CIRCLE_RADIUS) { 
1898           typeres = IntAna_Point;
1899           Alpha = (R1 + (dO1O2 - R2)) * 0.5;
1900         }
1901         else { 
1902           typeres = IntAna_Circle;
1903           dir1 = Dir;
1904           param1 = Beta;
1905         }          
1906         pt1.SetCoord( O1.X() + Alpha*Dir.X()
1907                      ,O1.Y() + Alpha*Dir.Y()
1908                      ,O1.Z() + Alpha*Dir.Z());
1909         
1910         nbint=1;
1911       }
1912     }
1913   }
1914 }
1915
1916 //=======================================================================
1917 //function : IntAna_QuadQuadGeo
1918 //purpose  : Plane - Torus
1919 //=======================================================================
1920 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& Pln,
1921                                        const gp_Torus& Tor,
1922                                        const Standard_Real Tol) 
1923 : done(Standard_False),
1924   nbint(0),
1925   typeres(IntAna_Empty),
1926   pt1(0,0,0),
1927   pt2(0,0,0),
1928   pt3(0,0,0),
1929   pt4(0,0,0),
1930   param1(0),
1931   param2(0),
1932   param3(0),
1933   param4(0),
1934   param1bis(0),
1935   param2bis(0),
1936   myCommonGen(Standard_False),
1937   myPChar(0,0,0)
1938 {
1939   InitTolerances();
1940   Perform(Pln,Tor,Tol);
1941 }
1942 //=======================================================================
1943 //function : Perform
1944 //purpose  : 
1945 //=======================================================================
1946 void IntAna_QuadQuadGeo::Perform(const gp_Pln& Pln,
1947                                  const gp_Torus& Tor,
1948                                  const Standard_Real Tol)
1949 {
1950   done = Standard_True;
1951   //
1952   Standard_Real aRMin, aRMaj;
1953   //
1954   aRMin = Tor.MinorRadius();
1955   aRMaj = Tor.MajorRadius();
1956   if (aRMin >= aRMaj) {
1957     typeres = IntAna_NoGeometricSolution; 
1958     return;
1959   }
1960   //
1961   const gp_Ax1 aPlnAx = Pln.Axis();
1962   const gp_Ax1 aTorAx = Tor.Axis();
1963   //
1964   Standard_Boolean bParallel, bNormal;
1965   //
1966   bParallel = aTorAx.IsParallel(aPlnAx, myEPSILON_AXES_PARA);
1967   bNormal = !bParallel ? aTorAx.IsNormal(aPlnAx, myEPSILON_AXES_PARA) : Standard_False;
1968   if (!bNormal && !bParallel) {
1969     typeres = IntAna_NoGeometricSolution; 
1970     return;
1971   }
1972   //
1973   Standard_Real aDist;
1974   //
1975   gp_Pnt aTorLoc = aTorAx.Location();
1976   if (bParallel) {
1977     Standard_Real aDt, X, Y, Z, A, B, C, D;
1978     //
1979     Pln.Coefficients(A,B,C,D);
1980     aTorLoc.Coord(X,Y,Z);
1981     aDist = A*X + B*Y + C*Z + D;
1982     //
1983     if ((Abs(aDist) - aRMin) > Tol) {
1984       typeres=IntAna_Empty;
1985       return;
1986     }
1987     //
1988     typeres = IntAna_Circle;
1989     //
1990     pt1.SetCoord(X - aDist*A, Y - aDist*B, Z - aDist*C);
1991     aDt = Sqrt(Abs(aRMin*aRMin - aDist*aDist));
1992     param1 = aRMaj + aDt;
1993     dir1 = aTorAx.Direction();
1994     nbint = 1;
1995     if ((Abs(aDist) < aRMin) && (aDt > Tol)) {
1996       pt2 = pt1;
1997       param2 = aRMaj - aDt;
1998       dir2 = dir1;
1999       nbint = 2;
2000     }
2001   }
2002   //
2003   else {
2004     aDist = Pln.Distance(aTorLoc);
2005     if (aDist > myEPSILON_DISTANCE) {
2006       typeres = IntAna_NoGeometricSolution; 
2007       return;
2008     }
2009     //
2010     typeres = IntAna_Circle;
2011     param2 = param1 = aRMin;
2012     dir2 = dir1 = aPlnAx.Direction();
2013     nbint = 2;
2014     //
2015     gp_Dir aDir = aTorAx.Direction()^dir1;
2016     pt1.SetXYZ(aTorLoc.XYZ() + aRMaj*aDir.XYZ());
2017     pt2.SetXYZ(aTorLoc.XYZ() - aRMaj*aDir.XYZ());
2018   }
2019 }
2020
2021 //=======================================================================
2022 //function : IntAna_QuadQuadGeo
2023 //purpose  : Cylinder - Torus
2024 //=======================================================================
2025 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
2026                                        const gp_Torus& Tor,
2027                                        const Standard_Real Tol) 
2028 : done(Standard_False),
2029   nbint(0),
2030   typeres(IntAna_Empty),
2031   pt1(0,0,0),
2032   pt2(0,0,0),
2033   pt3(0,0,0),
2034   pt4(0,0,0),
2035   param1(0),
2036   param2(0),
2037   param3(0),
2038   param4(0),
2039   param1bis(0),
2040   param2bis(0),
2041   myCommonGen(Standard_False),
2042   myPChar(0,0,0)
2043 {
2044   InitTolerances();
2045   Perform(Cyl,Tor,Tol);
2046 }
2047 //=======================================================================
2048 //function : Perform
2049 //purpose  : 
2050 //=======================================================================
2051 void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl,
2052                                  const gp_Torus& Tor,
2053                                  const Standard_Real Tol) 
2054 {
2055   done = Standard_True;
2056   //
2057   Standard_Real aRMin, aRMaj;
2058   //
2059   aRMin = Tor.MinorRadius();
2060   aRMaj = Tor.MajorRadius();
2061   if (aRMin >= aRMaj) {
2062     typeres = IntAna_NoGeometricSolution; 
2063     return;
2064   }
2065   //
2066   const gp_Ax1 aCylAx = Cyl.Axis();
2067   const gp_Ax1 aTorAx = Tor.Axis();
2068   //
2069   const gp_Lin aLin(aTorAx);
2070   const gp_Pnt aLocCyl = Cyl.Location();
2071   //
2072   if (!aTorAx.IsParallel(aCylAx, myEPSILON_AXES_PARA) ||
2073       (aLin.Distance(aLocCyl) > myEPSILON_DISTANCE)) {
2074     typeres = IntAna_NoGeometricSolution; 
2075     return;
2076   }
2077   //
2078   Standard_Real aRCyl;
2079   //
2080   aRCyl = Cyl.Radius();
2081   if (((aRCyl + Tol) < (aRMaj - aRMin)) || ((aRCyl - Tol) > (aRMaj + aRMin))) {
2082     typeres = IntAna_Empty;
2083     return;
2084   }
2085   //
2086   typeres = IntAna_Circle;
2087   //
2088   Standard_Real aDist = Sqrt(Abs(aRMin*aRMin - (aRCyl-aRMaj)*(aRCyl-aRMaj)));
2089   gp_XYZ aTorLoc = aTorAx.Location().XYZ();
2090   //
2091   dir1 = aTorAx.Direction();
2092   pt1.SetXYZ(aTorLoc + aDist*dir1.XYZ());
2093   param1 = aRCyl;
2094   nbint = 1;
2095   if ((aDist > Tol) && (aRCyl > (aRMaj - aRMin)) &&
2096                        (aRCyl < (aRMaj + aRMin))) {
2097     dir2 = dir1;
2098     pt2.SetXYZ(aTorLoc - aDist*dir2.XYZ());
2099     param2 = param1;
2100     nbint = 2;
2101   }
2102 }
2103
2104 //=======================================================================
2105 //function : IntAna_QuadQuadGeo
2106 //purpose  : Cone - Torus
2107 //=======================================================================
2108 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cone& Con,
2109                                        const gp_Torus& Tor,
2110                                        const Standard_Real Tol) 
2111 : done(Standard_False),
2112   nbint(0),
2113   typeres(IntAna_Empty),
2114   pt1(0,0,0),
2115   pt2(0,0,0),
2116   pt3(0,0,0),
2117   pt4(0,0,0),
2118   param1(0),
2119   param2(0),
2120   param3(0),
2121   param4(0),
2122   param1bis(0),
2123   param2bis(0),
2124   myCommonGen(Standard_False),
2125   myPChar(0,0,0)
2126 {
2127   InitTolerances();
2128   Perform(Con,Tor,Tol);
2129 }
2130 //=======================================================================
2131 //function : Perform
2132 //purpose  : 
2133 //=======================================================================
2134 void IntAna_QuadQuadGeo::Perform(const gp_Cone& Con,
2135                                  const gp_Torus& Tor,
2136                                  const Standard_Real Tol) 
2137 {
2138   done = Standard_True;
2139   //
2140   Standard_Real aRMin, aRMaj;
2141   //
2142   aRMin = Tor.MinorRadius();
2143   aRMaj = Tor.MajorRadius();
2144   if (aRMin >= aRMaj) {
2145     typeres = IntAna_NoGeometricSolution; 
2146     return;
2147   }
2148   //
2149   const gp_Ax1 aConAx = Con.Axis();
2150   const gp_Ax1 aTorAx = Tor.Axis();
2151   //
2152   const gp_Lin aLin(aTorAx);
2153   const gp_Pnt aConApex = Con.Apex();
2154   //
2155   if (!aTorAx.IsParallel(aConAx, myEPSILON_AXES_PARA) ||
2156       (aLin.Distance(aConApex) > myEPSILON_DISTANCE)) {
2157     typeres = IntAna_NoGeometricSolution; 
2158     return;
2159   }
2160   //
2161   Standard_Real anAngle, aDist, aParam[4], aDt;
2162   Standard_Integer i;
2163   gp_Pnt aTorLoc, aPCT, aPN, aPt[4];
2164   gp_Dir aDir[4];
2165   //
2166   anAngle = Con.SemiAngle();
2167   aTorLoc = aTorAx.Location();
2168   //
2169   aPN.SetXYZ(aTorLoc.XYZ() + aRMaj*Tor.YAxis().Direction().XYZ());
2170   gp_Dir aDN (gp_Vec(aTorLoc, aPN));
2171   gp_Ax1 anAxCLRot(aConApex, aDN);
2172   gp_Lin aConL = aLin.Rotated(anAxCLRot, anAngle);
2173   gp_Dir aDL = aConL.Position().Direction();
2174   gp_Dir aXDir = Tor.XAxis().Direction();
2175   //
2176   typeres = IntAna_Empty;
2177   //
2178   for (i = 0; i < 2; ++i) {
2179     if (i) {
2180       aXDir.Reverse();
2181     }
2182     aPCT.SetXYZ(aTorLoc.XYZ() + aRMaj*aXDir.XYZ());
2183     //
2184     aDist = aConL.Distance(aPCT);
2185     if (aDist > aRMin+Tol) {
2186       continue;
2187     }
2188     //
2189     typeres = IntAna_Circle;
2190     //
2191     gp_XYZ aPh = aPCT.XYZ() - aDist*aConL.Normal(aPCT).Direction().XYZ();
2192     aDt = Sqrt(Abs(aRMin*aRMin - aDist*aDist));
2193     //
2194     gp_Pnt aP;
2195     gp_XYZ aDVal = aDt*aDL.XYZ();
2196     aP.SetXYZ(aPh + aDVal);
2197     aParam[nbint] = aLin.Distance(aP);
2198     aPt[nbint].SetXYZ(aP.XYZ() - aParam[nbint]*aXDir.XYZ());
2199     aDir[nbint] = aTorAx.Direction();
2200     ++nbint;
2201     if ((aDist < aRMin) && (aDt > Tol)) {
2202       aP.SetXYZ(aPh - aDVal);
2203       aParam[nbint] = aLin.Distance(aP);
2204       aPt[nbint].SetXYZ(aP.XYZ() - aParam[nbint]*aXDir.XYZ());
2205       aDir[nbint] = aDir[nbint-1];
2206       ++nbint;
2207     }
2208   }
2209   //
2210   for (i = 0; i < nbint; ++i) {
2211     switch (i) {
2212     case 0:{
2213       pt1 = aPt[i];
2214       param1 = aParam[i];
2215       dir1 = aDir[i];
2216       break;
2217     }
2218     case 1:{
2219       pt2 = aPt[i];
2220       param2 = aParam[i];
2221       dir2 = aDir[i];
2222       break;
2223     }
2224     case 2:{
2225       pt3 = aPt[i];
2226       param3 = aParam[i];
2227       dir3 = aDir[i];
2228       break;
2229     }
2230     case 3:{
2231       pt4 = aPt[i];
2232       param4 = aParam[i];
2233       dir4 = aDir[i];
2234       break;
2235     }
2236     default:
2237       break;
2238     }
2239   }
2240 }
2241
2242 //=======================================================================
2243 //function : IntAna_QuadQuadGeo
2244 //purpose  : Sphere - Torus
2245 //=======================================================================
2246 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Sphere& Sph,
2247                                        const gp_Torus& Tor,
2248                                        const Standard_Real Tol) 
2249 : done(Standard_False),
2250   nbint(0),
2251   typeres(IntAna_Empty),
2252   pt1(0,0,0),
2253   pt2(0,0,0),
2254   pt3(0,0,0),
2255   pt4(0,0,0),
2256   param1(0),
2257   param2(0),
2258   param3(0),
2259   param4(0),
2260   param1bis(0),
2261   param2bis(0),
2262   myCommonGen(Standard_False),
2263   myPChar(0,0,0)
2264 {
2265   InitTolerances();
2266   Perform(Sph,Tor,Tol);
2267 }
2268 //=======================================================================
2269 //function : Perform
2270 //purpose  : 
2271 //=======================================================================
2272 void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph,
2273                                  const gp_Torus& Tor,
2274                                  const Standard_Real Tol) 
2275 {
2276   done = Standard_True;
2277   //
2278   Standard_Real aRMin, aRMaj;
2279   //
2280   aRMin = Tor.MinorRadius();
2281   aRMaj = Tor.MajorRadius();
2282   if (aRMin >= aRMaj) {
2283     typeres = IntAna_NoGeometricSolution; 
2284     return;
2285   }
2286   //
2287   const gp_Ax1 aTorAx = Tor.Axis();
2288   const gp_Lin aLin(aTorAx);
2289   const gp_Pnt aSphLoc = Sph.Location();
2290   //
2291   if (aLin.Distance(aSphLoc) > myEPSILON_DISTANCE) {
2292     typeres = IntAna_NoGeometricSolution;
2293     return;
2294   }
2295   //
2296   Standard_Real aRSph, aDist;
2297   gp_Pnt aTorLoc;
2298   //
2299   gp_Dir aXDir = Tor.XAxis().Direction();
2300   aTorLoc.SetXYZ(aTorAx.Location().XYZ() + aRMaj*aXDir.XYZ());
2301   aRSph = Sph.Radius();
2302   //
2303   gp_Vec aVec12(aTorLoc, aSphLoc);
2304   aDist = aVec12.Magnitude();
2305   if (((aDist - Tol) > (aRMin + aRSph)) || 
2306       ((aDist + Tol) < Abs(aRMin - aRSph))) {
2307     typeres = IntAna_Empty;
2308     return;
2309   }
2310   //
2311   typeres = IntAna_Circle;
2312   //
2313   Standard_Real anAlpha, aBeta;
2314   //
2315   anAlpha = 0.5*(aRMin*aRMin - aRSph*aRSph + aDist*aDist ) / aDist;
2316   aBeta = Sqrt(Abs(aRMin*aRMin - anAlpha*anAlpha));
2317   //
2318   gp_Dir aDir12(aVec12);
2319   gp_XYZ aPh = aTorLoc.XYZ() + anAlpha*aDir12.XYZ();
2320   gp_Dir aDC = Tor.YAxis().Direction()^aDir12;
2321   //
2322   gp_Pnt aP;
2323   gp_XYZ aDVal = aBeta*aDC.XYZ();
2324   aP.SetXYZ(aPh + aDVal);
2325   param1 = aLin.Distance(aP);
2326   pt1.SetXYZ(aP.XYZ() - param1*aXDir.XYZ());
2327   dir1 = aTorAx.Direction();
2328   nbint = 1;
2329   if ((aDist < (aRSph + aRMin)) && (aDist > Abs(aRSph - aRMin)) && 
2330       (aDVal.Modulus() > Tol)) {
2331     aP.SetXYZ(aPh - aDVal);
2332     param2 = aLin.Distance(aP);
2333     pt2.SetXYZ(aP.XYZ() - param2*aXDir.XYZ());
2334     dir2 = dir1;
2335     nbint = 2;
2336   }
2337 }
2338
2339 //=======================================================================
2340 //function : IntAna_QuadQuadGeo
2341 //purpose  : Torus - Torus
2342 //=======================================================================
2343 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Torus& Tor1,
2344                                        const gp_Torus& Tor2,
2345                                        const Standard_Real Tol) 
2346 : done(Standard_False),
2347   nbint(0),
2348   typeres(IntAna_Empty),
2349   pt1(0,0,0),
2350   pt2(0,0,0),
2351   pt3(0,0,0),
2352   pt4(0,0,0),
2353   param1(0),
2354   param2(0),
2355   param3(0),
2356   param4(0),
2357   param1bis(0),
2358   param2bis(0),
2359   myCommonGen(Standard_False),
2360   myPChar(0,0,0)
2361 {
2362   InitTolerances();
2363   Perform(Tor1,Tor2,Tol);
2364 }
2365 //=======================================================================
2366 //function : Perform
2367 //purpose  : 
2368 //=======================================================================
2369 void IntAna_QuadQuadGeo::Perform(const gp_Torus& Tor1,
2370                                  const gp_Torus& Tor2,
2371                                  const Standard_Real Tol) 
2372 {
2373   done = Standard_True;
2374   //
2375   Standard_Real aRMin1, aRMin2, aRMaj1, aRMaj2;
2376   //
2377   aRMin1 = Tor1.MinorRadius();
2378   aRMaj1 = Tor1.MajorRadius();
2379   aRMin2 = Tor2.MinorRadius();
2380   aRMaj2 = Tor2.MajorRadius();
2381   if (aRMin1 >= aRMaj1 || aRMin2 >= aRMaj2) {
2382     typeres = IntAna_NoGeometricSolution;
2383     return;
2384   }
2385   //
2386   const gp_Ax1 anAx1 = Tor1.Axis();
2387   const gp_Ax1 anAx2 = Tor2.Axis();
2388   //
2389   gp_Lin aL1(anAx1);
2390   if (!anAx1.IsParallel(anAx2, myEPSILON_AXES_PARA) ||
2391       (aL1.Distance(anAx2.Location()) > myEPSILON_DISTANCE)) {
2392     typeres = IntAna_NoGeometricSolution;
2393     return;
2394   }
2395   //
2396   gp_Pnt aLoc1, aLoc2;
2397   //
2398   aLoc1 = anAx1.Location();
2399   aLoc2 = anAx2.Location();
2400   //
2401   if (aLoc1.IsEqual(aLoc2, Tol) &&
2402       (Abs(aRMin1 - aRMin2) <= Tol) && 
2403       (Abs(aRMaj1 - aRMaj2) <= Tol)) {
2404     typeres = IntAna_Same;
2405     return;
2406   }
2407   //
2408   Standard_Real aDist;
2409   gp_Pnt aP1, aP2;
2410   //
2411   gp_Dir aXDir1 = Tor1.XAxis().Direction();
2412   aP1.SetXYZ(aLoc1.XYZ() + aRMaj1*aXDir1.XYZ());
2413   aP2.SetXYZ(aLoc2.XYZ() + aRMaj2*aXDir1.XYZ());
2414   //
2415   gp_Vec aV12(aP1, aP2);
2416   aDist = aV12.Magnitude();
2417   if (((aDist - Tol) > (aRMin1 + aRMin2)) || 
2418       ((aDist + Tol) < Abs(aRMin1 - aRMin2))) {
2419     typeres = IntAna_Empty;
2420     return;
2421   }
2422   //
2423   typeres = IntAna_Circle;
2424   //
2425   Standard_Real anAlpha, aBeta;
2426   //
2427   anAlpha = 0.5*(aRMin1*aRMin1 - aRMin2*aRMin2 + aDist*aDist ) / aDist;
2428   aBeta = Sqrt(Abs(aRMin1*aRMin1 - anAlpha*anAlpha));
2429   //
2430   gp_Dir aDir12(aV12);
2431   gp_XYZ aPh = aP1.XYZ() + anAlpha*aDir12.XYZ();
2432   gp_Dir aDC = Tor1.YAxis().Direction()^aDir12;
2433   //
2434   gp_Pnt aP;
2435   gp_XYZ aDVal = aBeta*aDC.XYZ();
2436   aP.SetXYZ(aPh + aDVal);
2437   param1 = aL1.Distance(aP);
2438   pt1.SetXYZ(aP.XYZ() - param1*aXDir1.XYZ());
2439   dir1 = anAx1.Direction();
2440   nbint = 1;
2441   if ((aDist < (aRMin1 + aRMin2)) && (aDist > Abs(aRMin1 - aRMin2)) && 
2442       aDVal.Modulus() > Tol) {
2443     aP.SetXYZ(aPh - aDVal);
2444     param2 = aL1.Distance(aP);
2445     pt2.SetXYZ(aP.XYZ() - param2*aXDir1.XYZ());
2446     dir2 = dir1;
2447     nbint = 2;
2448   }
2449 }
2450
2451 //=======================================================================
2452 //function : Point
2453 //purpose  : Returns a Point
2454 //=======================================================================
2455   gp_Pnt IntAna_QuadQuadGeo::Point(const Standard_Integer n) const 
2456 {
2457   if(!done)          {    StdFail_NotDone::Raise();        }
2458   if(n>nbint || n<1) {    Standard_DomainError::Raise();   }
2459   if(typeres==IntAna_PointAndCircle) {
2460     if(n!=1) { Standard_DomainError::Raise();  }
2461     if(param1==0.0) return(pt1);
2462     return(pt2);
2463   }
2464   else if(typeres==IntAna_Point) {
2465     if(n==1) return(pt1);
2466     return(pt2);
2467   }
2468
2469   // WNT (what can you expect from MicroSoft ?)
2470   return gp_Pnt(0,0,0);
2471 }
2472 //=======================================================================
2473 //function : Line
2474 //purpose  : Returns a Line
2475 //=======================================================================
2476   gp_Lin IntAna_QuadQuadGeo::Line(const Standard_Integer n) const 
2477 {
2478   if(!done)        {   StdFail_NotDone::Raise();   }
2479   if((n>nbint) || (n<1) || (typeres!=IntAna_Line)) {
2480     Standard_DomainError::Raise();
2481     }
2482   if(n==1) {  return(gp_Lin(pt1,dir1));   }
2483   else {      return(gp_Lin(pt2,dir2));   }
2484 }
2485 //=======================================================================
2486 //function : Circle
2487 //purpose  : Returns a Circle
2488 //=======================================================================
2489   gp_Circ IntAna_QuadQuadGeo::Circle(const Standard_Integer n) const 
2490 {
2491   if(!done) {    StdFail_NotDone::Raise();     }
2492   if(typeres==IntAna_PointAndCircle) {
2493     if(n!=1) { Standard_DomainError::Raise();  }
2494     if(param2==0.0) return(gp_Circ(DirToAx2(pt1,dir1),param1));
2495     return(gp_Circ(DirToAx2(pt2,dir2),param2));
2496   }
2497   else if((n>nbint) || (n<1) || (typeres!=IntAna_Circle)) {
2498     Standard_DomainError::Raise();
2499     }
2500   if      (n==1) { return(gp_Circ(DirToAx2(pt1,dir1),param1));}
2501   else if (n==2) { return(gp_Circ(DirToAx2(pt2,dir2),param2));}
2502   else if (n==3) { return(gp_Circ(DirToAx2(pt3,dir3),param3));}
2503   else           { return(gp_Circ(DirToAx2(pt4,dir4),param4));}
2504 }
2505
2506 //=======================================================================
2507 //function : Ellipse
2508 //purpose  : Returns a Elips  
2509 //=======================================================================
2510   gp_Elips IntAna_QuadQuadGeo::Ellipse(const Standard_Integer n) const
2511 {
2512   if(!done) {     StdFail_NotDone::Raise();     }
2513   if((n>nbint) || (n<1) || (typeres!=IntAna_Ellipse)) {
2514     Standard_DomainError::Raise();
2515   }
2516
2517   if(n==1) {
2518     Standard_Real R1=param1, R2=param1bis, aTmp;
2519     if (R1<R2) {
2520       aTmp=R1; R1=R2; R2=aTmp;
2521     }
2522     gp_Ax2 anAx2(pt1, dir1 ,dir2);
2523     gp_Elips anElips (anAx2, R1, R2);
2524     return anElips;
2525   }
2526   else {
2527     Standard_Real R1=param2, R2=param2bis, aTmp;
2528     if (R1<R2) {
2529       aTmp=R1; R1=R2; R2=aTmp;
2530     }
2531     gp_Ax2 anAx2(pt2, dir2 ,dir1);
2532     gp_Elips anElips (anAx2, R1, R2);
2533     return anElips;
2534   }
2535 }
2536 //=======================================================================
2537 //function : Parabola
2538 //purpose  : Returns a Parabola 
2539 //=======================================================================
2540   gp_Parab IntAna_QuadQuadGeo::Parabola(const Standard_Integer n) const 
2541 {
2542   if(!done) {
2543     StdFail_NotDone::Raise();
2544     }
2545   if (typeres!=IntAna_Parabola) {
2546     Standard_DomainError::Raise();
2547   }
2548   if((n>nbint) || (n!=1)) {
2549     Standard_OutOfRange::Raise();
2550   }
2551   return(gp_Parab(gp_Ax2( pt1
2552                          ,dir1
2553                          ,dir2)
2554                   ,param1));
2555 }
2556 //=======================================================================
2557 //function : Hyperbola
2558 //purpose  : Returns a Hyperbola  
2559 //=======================================================================
2560   gp_Hypr IntAna_QuadQuadGeo::Hyperbola(const Standard_Integer n) const 
2561 {
2562   if(!done) {
2563     StdFail_NotDone::Raise();
2564     }
2565   if((n>nbint) || (n<1) || (typeres!=IntAna_Hyperbola)) {
2566     Standard_DomainError::Raise();
2567     }
2568   if(n==1) {
2569     return(gp_Hypr(gp_Ax2( pt1
2570                           ,dir1
2571                           ,dir2)
2572                    ,param1,param1bis));
2573   }
2574   else {
2575     return(gp_Hypr(gp_Ax2( pt2
2576                           ,dir1
2577                           ,dir2.Reversed())
2578                    ,param2,param2bis));
2579   }
2580 }
2581 //=======================================================================
2582 //function : HasCommonGen
2583 //purpose  : 
2584 //=======================================================================
2585 Standard_Boolean IntAna_QuadQuadGeo::HasCommonGen() const
2586 {
2587   return myCommonGen;
2588 }
2589 //=======================================================================
2590 //function : PChar
2591 //purpose  : 
2592 //=======================================================================
2593 const gp_Pnt& IntAna_QuadQuadGeo::PChar() const
2594 {
2595   return myPChar;
2596 }
2597 //=======================================================================
2598 //function : RefineDir
2599 //purpose  : 
2600 //=======================================================================
2601 void RefineDir(gp_Dir& aDir)
2602 {
2603   Standard_Integer k, m, n;
2604   Standard_Real aC[3];
2605   //
2606   aDir.Coord(aC[0], aC[1], aC[2]);
2607   //
2608   m=0;
2609   n=0;
2610   for (k=0; k<3; ++k) {
2611     if (aC[k]==1. || aC[k]==-1.) {
2612       ++m;
2613     }
2614     else if (aC[k]!=0.) {
2615       ++n;
2616     }
2617   }
2618   //
2619   if (m && n) {
2620     Standard_Real aEps, aR1, aR2, aNum;
2621     //
2622     aEps=RealEpsilon();
2623     aR1=1.-aEps;
2624     aR2=1.+aEps;
2625     //
2626     for (k=0; k<3; ++k) {
2627       m=(aC[k]>0.);
2628       aNum=(m)? aC[k] : -aC[k];
2629       if (aNum>aR1 && aNum<aR2) {
2630         if (m) {
2631           aC[k]=1.;
2632         }          
2633         else {
2634           aC[k]=-1.;
2635         }
2636         //
2637         aC[(k+1)%3]=0.;
2638         aC[(k+2)%3]=0.;
2639         break;
2640       }
2641     }
2642     aDir.SetCoord(aC[0], aC[1], aC[2]);
2643   }
2644 }
2645
2646
2647