e803a1ccb9333d27aa37d66a42789f468301d1ff
[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-2012 OPEN CASCADE SAS
5 //
6 // The content of this file is subject to the Open CASCADE Technology Public
7 // License Version 6.5 (the "License"). You may not use the content of this file
8 // except in compliance with the License. Please obtain a copy of the License
9 // at http://www.opencascade.org and read it completely before using this file.
10 //
11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
13 //
14 // The Original Code and all software distributed under the License is
15 // distributed on an "AS IS" basis, without warranty of any kind, and the
16 // Initial Developer hereby disclaims all such warranties, including without
17 // limitation, any warranties of merchantability, fitness for a particular
18 // purpose or non-infringement. Please see the License for the specific terms
19 // and conditions governing the rights and limitations under the License.
20
21 //----------------------------------------------------------------------
22 //-- Purposse: Geometric Intersection between two Natural Quadric 
23 //--          If the intersection is not a conic, 
24 //--          analytical methods must be called.
25 //----------------------------------------------------------------------
26 #ifndef DEB
27 #define No_Standard_RangeError
28 #define No_Standard_OutOfRange
29 #endif
30
31 #include <IntAna_QuadQuadGeo.ixx>
32
33 #include <IntAna_IntConicQuad.hxx>
34 #include <StdFail_NotDone.hxx>
35 #include <Standard_DomainError.hxx>
36 #include <Standard_OutOfRange.hxx>
37 #include <math_DirectPolynomialRoots.hxx>
38
39 #include <gp.hxx>
40 #include <gp_Pln.hxx>
41 #include <gp_Vec.hxx>
42 #include <ElSLib.hxx>
43 #include <ElCLib.hxx>
44
45 #include <gp_Dir.hxx>
46 #include <gp_XYZ.hxx>
47 #include <gp_Pnt2d.hxx>
48 #include <gp_Vec2d.hxx>
49 #include <gp_Dir2d.hxx>
50
51
52 static
53   gp_Ax2 DirToAx2(const gp_Pnt& P,const gp_Dir& D);
54 static
55   void RefineDir(gp_Dir& aDir);
56
57 //=======================================================================
58 //class :  
59 //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 
60 //=======================================================================
61 class AxeOperator {
62  public:
63   AxeOperator(const gp_Ax1& A1,const gp_Ax1& A2);
64
65   void Distance(Standard_Real& dist,
66                 Standard_Real& Param1,
67                 Standard_Real& Param2);
68   
69   gp_Pnt PtIntersect()              { 
70     return ptintersect;
71   }
72   Standard_Boolean Coplanar(void)   { 
73     return thecoplanar;
74   }
75   Standard_Boolean Same(void)       {
76     return theparallel && (thedistance<myEPSILON_DISTANCE); 
77   }
78   Standard_Real Distance(void)      { 
79     return thedistance ;
80   }
81   Standard_Boolean Intersect(void)  { 
82     return (thecoplanar && (!theparallel));
83   }
84   Standard_Boolean Parallel(void)   { 
85     return theparallel; 
86   }
87   Standard_Boolean Normal(void)     { 
88     return thenormal;
89   }
90
91  protected:
92   Standard_Real Det33(const Standard_Real a11,
93                       const Standard_Real a12,
94                       const Standard_Real a13,
95                       const Standard_Real a21,
96                       const Standard_Real a22,
97                       const Standard_Real a23,
98                       const Standard_Real a31,
99                       const Standard_Real a32,
100                       const Standard_Real a33) {
101     Standard_Real theReturn =  
102       a11*(a22*a33-a32*a23) - a21*(a12*a33-a32*a13) + a31*(a12*a23-a22*a13) ;   
103     return theReturn ;
104   }
105
106  private:
107   gp_Pnt ptintersect;
108   gp_Ax1 Axe1;
109   gp_Ax1 Axe2;
110   Standard_Real thedistance;
111   Standard_Boolean theparallel;
112   Standard_Boolean thecoplanar;
113   Standard_Boolean thenormal;
114   //
115   Standard_Real myEPSILON_DISTANCE;
116   Standard_Real myEPSILON_AXES_PARA;
117 };
118
119 //=======================================================================
120 //function : AxeOperator::AxeOperator
121 //purpose  : 
122 //=======================================================================
123   AxeOperator::AxeOperator(const gp_Ax1& A1,const gp_Ax1& A2) 
124 {
125   myEPSILON_DISTANCE=0.00000000000001;
126   myEPSILON_AXES_PARA=0.000000000001;
127   Axe1=A1; 
128   Axe2=A2;
129   //---------------------------------------------------------------------
130   gp_Dir V1=Axe1.Direction();
131   gp_Dir V2=Axe2.Direction();
132   gp_Pnt P1=Axe1.Location();
133   gp_Pnt P2=Axe2.Location();
134   //
135   RefineDir(V1);
136   RefineDir(V2);
137   thecoplanar= Standard_False;
138   thenormal  = Standard_False;
139
140   //--- check if the two axis are parallel
141   theparallel=V1.IsParallel(V2, myEPSILON_AXES_PARA);  
142   //--- Distance between the two axis
143   gp_XYZ perp(A1.Direction().XYZ().Crossed(A2.Direction().XYZ()));
144   if (theparallel) { 
145     gp_Lin L1(A1);
146     thedistance = L1.Distance(A2.Location());
147   }
148   else {   
149     thedistance = Abs(gp_Vec(perp.Normalized()).Dot(gp_Vec(Axe1.Location(),
150                                                            Axe2.Location())));
151   }
152   //--- check if Axis are Coplanar
153   Standard_Real D33;
154   if(thedistance<myEPSILON_DISTANCE) {
155     D33=Det33(V1.X(),V1.Y(),V1.Z()
156               ,V2.X(),V2.Y(),V2.Z()
157               ,P1.X()-P2.X(),P1.Y()-P2.Y(),P1.Z()-P2.Z());
158     if(Abs(D33)<=myEPSILON_DISTANCE) { 
159       thecoplanar=Standard_True;
160     }
161   }
162   else {
163     thecoplanar=Standard_True;
164     thenormal=(V1.Dot(V2)==0.0)? Standard_True : Standard_False;
165   }
166   //--- check if the two axis are concurrent
167   if(thecoplanar && (!theparallel)) {
168     Standard_Real smx=P2.X() - P1.X();
169     Standard_Real smy=P2.Y() - P1.Y();
170     Standard_Real smz=P2.Z() - P1.Z();
171     Standard_Real Det1,Det2,Det3,A;
172     Det1=V1.Y() * V2.X() - V1.X() * V2.Y();
173     Det2=V1.Z() * V2.Y() - V1.Y() * V2.Z();
174     Det3=V1.Z() * V2.X() - V1.X() * V2.Z();
175     
176     if((Det1!=0.0) && ((Abs(Det1) >= Abs(Det2))&&(Abs(Det1) >= Abs(Det3)))) {
177       A=(smy * V2.X() - smx * V2.Y())/Det1;
178     }
179     else if((Det2!=0.0) 
180             && ((Abs(Det2) >= Abs(Det1))
181                 &&(Abs(Det2) >= Abs(Det3)))) {
182       A=(smz * V2.Y() - smy * V2.Z())/Det2;
183     }
184     else {
185       A=(smz * V2.X() - smx * V2.Z())/Det3;
186     }
187     ptintersect.SetCoord( P1.X() + A * V1.X()
188                          ,P1.Y() + A * V1.Y()
189                          ,P1.Z() + A * V1.Z());
190   }
191   else { 
192     ptintersect.SetCoord(0,0,0);  //-- Pour eviter des FPE
193   }
194 }
195 //=======================================================================
196 //function : Distance
197 //purpose  : 
198 //=======================================================================
199   void AxeOperator::Distance(Standard_Real& dist,Standard_Real& Param1,Standard_Real& Param2)
200  {
201   gp_Vec O1O2(Axe1.Location(),Axe2.Location());   //-----------------------------
202   gp_Dir U1 = Axe1.Direction();   //-- juste pour voir. 
203   gp_Dir U2 = Axe2.Direction();
204   
205   gp_Dir N  = U1.Crossed(U2);
206   Standard_Real D = Det33(U1.X(),U2.X(),N.X(),
207                           U1.Y(),U2.Y(),N.Y(),
208                           U1.Z(),U2.Z(),N.Z());
209   if(D) { 
210     dist = Det33(U1.X(),U2.X(),O1O2.X(),
211                  U1.Y(),U2.Y(),O1O2.Y(),
212                  U1.Z(),U2.Z(),O1O2.Z()) / D;
213     Param1 = Det33(O1O2.X(),U2.X(),N.X(),
214                    O1O2.Y(),U2.Y(),N.Y(),
215                    O1O2.Z(),U2.Z(),N.Z()) / (-D);
216     //------------------------------------------------------------
217     //-- On resout P1 * Dir1 + P2 * Dir2 + d * N = O1O2
218     //-- soit : Segment perpendiculaire : O1+P1 D1
219     //--                                  O2-P2 D2
220     Param2 = Det33(U1.X(),O1O2.X(),N.X(),
221                    U1.Y(),O1O2.Y(),N.Y(),
222                    U1.Z(),O1O2.Z(),N.Z()) / (D);
223   }
224 }
225 //=======================================================================
226 //function : DirToAx2
227 //purpose  : returns a gp_Ax2 where D is the main direction
228 //=======================================================================
229 gp_Ax2 DirToAx2(const gp_Pnt& P,const gp_Dir& D) 
230 {
231   Standard_Real x=D.X(); Standard_Real ax=Abs(x);
232   Standard_Real y=D.Y(); Standard_Real ay=Abs(y);
233   Standard_Real z=D.Z(); Standard_Real az=Abs(z);
234   if( (ax==0.0) || ((ax<ay) && (ax<az)) ) {
235     return(gp_Ax2(P,D,gp_Dir(gp_Vec(0.0,-z,y))));
236   }
237   else if( (ay==0.0) || ((ay<ax) && (ay<az)) ) {
238     return(gp_Ax2(P,D,gp_Dir(gp_Vec(-z,0.0,x))));
239   }
240   else {
241     return(gp_Ax2(P,D,gp_Dir(gp_Vec(-y,x,0.0))));
242   }
243 }
244 //=======================================================================
245 //function : IntAna_QuadQuadGeo
246 //purpose  : Empty constructor
247 //=======================================================================
248   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(void)
249     : done(Standard_False),
250       nbint(0),
251       typeres(IntAna_Empty),
252       pt1(0,0,0),
253       pt2(0,0,0),
254       param1(0),
255       param2(0),
256       param1bis(0),
257       param2bis(0),
258       myCommonGen(Standard_False),
259       myPChar(0,0,0)
260 {
261   InitTolerances();
262 }
263 //=======================================================================
264 //function : InitTolerances
265 //purpose  : 
266 //=======================================================================
267   void IntAna_QuadQuadGeo::InitTolerances()
268 {
269   myEPSILON_DISTANCE               = 0.00000000000001;
270   myEPSILON_ANGLE_CONE             = 0.000000000001;
271   myEPSILON_MINI_CIRCLE_RADIUS     = 0.000000001;
272   myEPSILON_CYLINDER_DELTA_RADIUS  = 0.0000000000001;
273   myEPSILON_CYLINDER_DELTA_DISTANCE= 0.0000001;
274   myEPSILON_AXES_PARA              = 0.000000000001;
275 }
276 //=======================================================================
277 //function : IntAna_QuadQuadGeo
278 //purpose  : Pln  Pln 
279 //=======================================================================
280   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P1, 
281                                          const gp_Pln& P2,
282                                          const Standard_Real TolAng,
283                                          const Standard_Real Tol)
284 : done(Standard_False),
285   nbint(0),
286   typeres(IntAna_Empty),
287   pt1(0,0,0),
288   pt2(0,0,0),
289   param1(0),
290   param2(0),
291   param1bis(0),
292   param2bis(0),
293   myCommonGen(Standard_False),
294   myPChar(0,0,0)
295 {
296   InitTolerances();
297   Perform(P1,P2,TolAng,Tol);
298 }
299 //=======================================================================
300 //function : Perform
301 //purpose  : 
302 //=======================================================================
303   void IntAna_QuadQuadGeo::Perform (const gp_Pln& P1, 
304                                     const gp_Pln& P2,
305                                     const Standard_Real TolAng,
306                                     const Standard_Real Tol)
307 {
308   done=Standard_False;
309   //
310   param2bis=0.0;
311
312   Standard_Real A1 = 0., B1 = 0., C1 = 0., D1 = 0., A2 = 0., B2 = 0., C2 = 0., D2 = 0.;
313   P1.Coefficients(A1,B1,C1,D1);
314   P2.Coefficients(A2,B2,C2,D2);
315   
316   gp_Vec vd(gp_Vec(A1,B1,C1).Crossed(gp_Vec(A2,B2,C2)));
317   Standard_Real dist1= A2*P1.Location().X() + B2*P1.Location().Y() + C2*P1.Location().Z() + D2;
318   Standard_Real dist2= A1*P2.Location().X() + B1*P2.Location().Y() + C1*P2.Location().Z() + D1;
319
320   if(vd.Magnitude() <=TolAng) {
321     // normalles are collinear - planes are same or parallel
322     typeres = (Abs(dist1) <= Tol && Abs(dist2) <= Tol) ? IntAna_Same : IntAna_Empty;
323   }
324   else {
325     Standard_Real denom=A1*A2 + B1*B2 + C1*C2;
326
327     Standard_Real denom2 = denom*denom;
328     Standard_Real ddenom = 1. - denom2;
329     //denom = ( Abs(ddenom) <= 1.e-9 ) ? 1.e-9 : ddenom;
330     denom = ( Abs(ddenom) <= 1.e-16 ) ? 1.e-16 : ddenom;
331       
332     Standard_Real par1 = dist1/denom;
333     Standard_Real par2 = -dist2/denom;
334       
335     gp_Vec inter1(gp_Vec(A1,B1,C1).Crossed(vd));
336     gp_Vec inter2(gp_Vec(A2,B2,C2).Crossed(vd));
337       
338     Standard_Real X1=P1.Location().X() + par1*inter1.X();
339     Standard_Real Y1=P1.Location().Y() + par1*inter1.Y();
340     Standard_Real Z1=P1.Location().Z() + par1*inter1.Z();
341     Standard_Real X2=P2.Location().X() + par2*inter2.X();
342     Standard_Real Y2=P2.Location().Y() + par2*inter2.Y();
343     Standard_Real Z2=P2.Location().Z() + par2*inter2.Z();
344       
345     pt1=gp_Pnt((X1+X2)*0.5, (Y1+Y2)*0.5, (Z1+Z2)*0.5);
346     dir1 = gp_Dir(vd);
347     typeres = IntAna_Line;
348     nbint = 1;
349  
350   }
351   done=Standard_True;
352 }
353 //=======================================================================
354 //function : IntAna_QuadQuadGeo
355 //purpose  : Pln Cylinder
356 //=======================================================================
357   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo( const gp_Pln& P
358        ,const gp_Cylinder& Cl
359        ,const Standard_Real Tolang
360        ,const Standard_Real Tol
361        ,const Standard_Real H)
362     : done(Standard_False),
363       nbint(0),
364       typeres(IntAna_Empty),
365       pt1(0,0,0),
366       pt2(0,0,0),
367       param1(0),
368       param2(0),
369       param1bis(0),
370       param2bis(0),
371       myCommonGen(Standard_False),
372       myPChar(0,0,0)
373 {
374   InitTolerances();
375   Perform(P,Cl,Tolang,Tol,H);
376 }
377 //=======================================================================
378 //function : Perform
379 //purpose  : 
380 //=======================================================================
381   void IntAna_QuadQuadGeo::Perform( const gp_Pln& P
382                                    ,const gp_Cylinder& Cl
383                                    ,const Standard_Real Tolang
384                                    ,const Standard_Real Tol
385                                    ,const Standard_Real H) 
386 {
387   done = Standard_False;
388   Standard_Real dist,radius;
389   Standard_Real A,B,C,D;
390   Standard_Real X,Y,Z;
391   Standard_Real sint,cost,h;
392   gp_XYZ axex,axey,omega;
393
394   
395   param2bis=0.0;
396   radius = Cl.Radius();
397
398   gp_Lin axec(Cl.Axis());
399   gp_XYZ normp(P.Axis().Direction().XYZ());
400
401   P.Coefficients(A,B,C,D);
402   axec.Location().Coord(X,Y,Z);
403   dist = A*X + B*Y + C*Z + D; // la distance axe/plan est evaluee a l origine de l axe.
404
405   Standard_Real tolang = Tolang;
406   Standard_Boolean newparams = Standard_False;
407
408   gp_Vec ldv( axec.Direction() );
409   gp_Vec npv( normp );
410   Standard_Real dA = Abs( ldv.Angle( npv ) );
411   if( dA > (M_PI/4.) )
412     {
413       Standard_Real dang = Abs( ldv.Angle( npv ) ) - M_PI/2.;
414       Standard_Real dangle = Abs( dang );
415       if( dangle > Tolang )
416         {
417           Standard_Real sinda = Abs( Sin( dangle ) );
418           Standard_Real dif = Abs( sinda - Tol );
419           if( dif < Tol )
420             {
421               tolang = sinda * 2.;
422               newparams = Standard_True;
423             }
424         }
425     }
426
427   nbint = 0;
428   IntAna_IntConicQuad inter(axec,P,tolang,Tol,H);
429
430   if (inter.IsParallel()) {
431     // Le resultat de l intersection Plan-Cylindre est de type droite.
432     // il y a 1 ou 2 droites
433
434     typeres = IntAna_Line;
435     omega.SetCoord(X-dist*A,Y-dist*B,Z-dist*C);
436
437     if (Abs(Abs(dist)-radius) < Tol)
438       {
439         nbint = 1;
440         pt1.SetXYZ(omega);
441
442         if( newparams )
443           {
444             gp_XYZ omegaXYZ(X,Y,Z);
445             gp_XYZ omegaXYZtrnsl( omegaXYZ + 100.*axec.Direction().XYZ() );
446             Standard_Real Xt,Yt,Zt,distt;
447             omegaXYZtrnsl.Coord(Xt,Yt,Zt);
448             distt = A*Xt + B*Yt + C*Zt + D;
449             gp_XYZ omega1( omegaXYZtrnsl.X()-distt*A, omegaXYZtrnsl.Y()-distt*B, omegaXYZtrnsl.Z()-distt*C );
450             gp_Pnt ppt1;
451             ppt1.SetXYZ( omega1 );
452             gp_Vec vv1(pt1,ppt1);
453             gp_Dir dd1( vv1 );
454             dir1 = dd1;
455           }
456         else
457           dir1 = axec.Direction();
458     }
459     else if (Abs(dist) < radius)
460       {
461         nbint = 2;
462         h = Sqrt(radius*radius - dist*dist);
463         axey = axec.Direction().XYZ().Crossed(normp); // axey est normalise
464
465         pt1.SetXYZ(omega - h*axey);
466         pt2.SetXYZ(omega + h*axey);
467
468         if( newparams )
469           { 
470             gp_XYZ omegaXYZ(X,Y,Z);
471             gp_XYZ omegaXYZtrnsl( omegaXYZ + 100.*axec.Direction().XYZ() );
472             Standard_Real Xt,Yt,Zt,distt,ht;
473             omegaXYZtrnsl.Coord(Xt,Yt,Zt);
474             distt = A*Xt + B*Yt + C*Zt + D;
475             //      ht = Sqrt(radius*radius - distt*distt);
476             Standard_Real anSqrtArg = radius*radius - distt*distt;
477             ht = (anSqrtArg > 0.) ? Sqrt(anSqrtArg) : 0.;
478
479             gp_XYZ omega1( omegaXYZtrnsl.X()-distt*A, omegaXYZtrnsl.Y()-distt*B, omegaXYZtrnsl.Z()-distt*C );
480             gp_Pnt ppt1,ppt2;
481             ppt1.SetXYZ( omega1 - ht*axey);
482             ppt2.SetXYZ( omega1 + ht*axey);
483             gp_Vec vv1(pt1,ppt1);
484             gp_Vec vv2(pt2,ppt2);
485             gp_Dir dd1( vv1 );
486             gp_Dir dd2( vv2 );
487             dir1 = dd1;
488             dir2 = dd2;
489           }
490         else
491           {
492             dir1 = axec.Direction();
493             dir2 = axec.Direction();
494           }
495     }
496     //  else nbint = 0
497
498     // debug JAG : le nbint = 0 doit etre remplace par typeres = IntAna_Empty
499     // et ne pas etre seulement supprime...
500
501     else {
502       typeres = IntAna_Empty;
503     }
504   }
505   else {     // Il y a un point d intersection. C est le centre du cercle
506              // ou de l ellipse solution.
507
508     nbint = 1;
509     axey = normp.Crossed(axec.Direction().XYZ());
510     sint = axey.Modulus();
511
512     pt1 = inter.Point(1);
513     
514     if (sint < Tol/radius) {
515
516       // on construit un cercle avec comme axes X et Y ceux du cylindre
517       typeres = IntAna_Circle;
518
519       dir1 = axec.Direction(); // axe Z
520       dir2 = Cl.Position().XDirection();
521       param1 = radius;
522     }
523     else {
524
525       // on construit un ellipse
526       typeres = IntAna_Ellipse;
527       cost = Abs(axec.Direction().XYZ().Dot(normp));
528       axex = axey.Crossed(normp);
529
530       dir1.SetXYZ(normp);         //Modif ds ce bloc 
531       dir2.SetXYZ(axex);
532
533       param1    = radius/cost;
534       param1bis = radius;
535     }
536   }
537   if(typeres == IntAna_Ellipse) { 
538     if(   param1>100000.0*param1bis 
539        || param1bis>100000.0*param1) { 
540       done = Standard_False;
541       return;
542     }
543   }
544   done = Standard_True;
545 }
546 //=======================================================================
547 //function : IntAna_QuadQuadGeo
548 //purpose  : Pln Cone
549 //=======================================================================
550   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P,
551                                          const gp_Cone& Co,
552                                          const Standard_Real Tolang,
553                                          const Standard_Real Tol) 
554 : done(Standard_False),
555   nbint(0),
556   typeres(IntAna_Empty),
557   pt1(0,0,0),
558   pt2(0,0,0),
559   param1(0),
560   param2(0),
561   param1bis(0),
562   param2bis(0),
563   myCommonGen(Standard_False),
564   myPChar(0,0,0)
565 {
566   InitTolerances();
567   Perform(P,Co,Tolang,Tol);
568 }
569 //=======================================================================
570 //function : Perform
571 //purpose  : 
572 //=======================================================================
573   void IntAna_QuadQuadGeo::Perform(const gp_Pln& P,
574                                    const gp_Cone& Co,
575                                    const Standard_Real Tolang,
576                                    const Standard_Real Tol) 
577 {
578
579   done = Standard_False;
580   nbint = 0;
581
582   Standard_Real A,B,C,D;
583   Standard_Real X,Y,Z;
584   Standard_Real dist,sint,cost,sina,cosa,angl,costa;
585   Standard_Real dh;
586   gp_XYZ axex,axey;
587
588   gp_Lin axec(Co.Axis());
589   P.Coefficients(A,B,C,D);
590   gp_Pnt apex(Co.Apex());
591
592   apex.Coord(X,Y,Z);
593   dist = A*X + B*Y + C*Z + D; // distance signee sommet du cone/ Plan
594
595   gp_XYZ normp = P.Axis().Direction().XYZ();
596   if(P.Direct()==Standard_False) {  //-- lbr le 14 jan 97
597     normp.Reverse();
598   }
599
600   axey = normp.Crossed(Co.Axis().Direction().XYZ());
601   axex = axey.Crossed(normp);
602
603
604   angl = Co.SemiAngle();
605
606   cosa = Cos(angl);
607   sina = Abs(Sin(angl));
608
609
610   // Angle entre la normale au plan et l axe du cone, ramene entre 0. et PI/2.
611
612   sint = axey.Modulus();
613   cost = Abs(Co.Axis().Direction().XYZ().Dot(normp));
614
615   // Le calcul de costa permet de determiner si le plan contient
616   // un generatrice du cone : on calcul Sin((PI/2. - t) - angl)
617
618   costa = cost*cosa - sint*sina;  // sin((PI/2 -t)-angl)=cos(t+angl)
619                                   // avec  t ramene entre 0 et pi/2.
620
621   if (Abs(dist) < Tol) {
622     // on considere que le plan contient le sommet du cone.
623     // les solutions possibles sont donc : 1 point, 1 droite, 2 droites
624     // selon l inclinaison du plan.
625
626     if (Abs(costa) < Tolang) {          // plan parallele a la generatrice
627       typeres = IntAna_Line;
628       nbint = 1;
629       gp_XYZ ptonaxe(apex.XYZ() + 10.*(Co.Axis().Direction().XYZ()));
630       // point sur l axe du cone cote z positif
631
632       dist = A*ptonaxe.X() + B*ptonaxe.Y() + C*ptonaxe.Z() + D;
633       ptonaxe = ptonaxe - dist*normp;
634       pt1 = apex;
635       dir1.SetXYZ(ptonaxe - pt1.XYZ());
636     }
637     else if (cost < sina) {   // plan "interieur" au cone
638       typeres = IntAna_Line;
639       nbint = 2;
640       pt1 = apex;
641       pt2 = apex;
642       dh = Sqrt(sina*sina-cost*cost)/cosa;
643       dir1.SetXYZ(axex + dh*axey);
644       dir2.SetXYZ(axex - dh*axey);
645     }
646     else {                              // plan "exterieur" au cone
647       typeres = IntAna_Point;
648       nbint = 1;
649       pt1 = apex;
650     }
651   }
652   else {
653     // Solutions possibles : cercle, ellipse, parabole, hyperbole selon
654     // l inclinaison du plan.
655     Standard_Real deltacenter, distance;
656
657     if (cost < Tolang) {
658       // Le plan contient la direction de l axe du cone. La solution est
659       // l hyperbole
660       typeres = IntAna_Hyperbola;
661       nbint = 2;
662       pt1.SetXYZ(apex.XYZ()-dist*normp);
663       pt2 = pt1;
664       dir1=normp;
665       dir2.SetXYZ(axex);
666       param1     = param2 = Abs(dist/Tan(angl));
667       param1bis  = param2bis = Abs(dist);
668     }
669     else {
670
671       IntAna_IntConicQuad inter(axec,P,Tolang); // on a necessairement 1 point.
672       
673       gp_Pnt center(inter.Point(1));
674
675       // En fonction de la position de l intersection par rapport au sommet
676       // du cone, on change l axe x en -x et y en -y. Le parametre du sommet
677       // sur axec est negatif (voir definition du cone)
678
679       distance = apex.Distance(center);
680
681       if (inter.ParamOnConic(1) + Co.RefRadius()/Tan(angl) < 0.) {
682         axex.Reverse();
683         axey.Reverse();
684       }
685
686       if (Abs(costa) < Tolang) {  // plan parallele a une generatrice
687         typeres = IntAna_Parabola;
688         nbint = 1;
689         deltacenter = distance/2./cosa;
690         axex.Normalize();
691         pt1.SetXYZ(center.XYZ()-deltacenter*axex);
692         dir1 = normp;
693         dir2.SetXYZ(axex);
694         param1 = deltacenter*sina*sina;
695       }
696       else if (sint  < Tolang) {            // plan perpendiculaire a l axe
697         typeres = IntAna_Circle;
698         nbint = 1;
699         pt1 = center;
700         dir1 = Co.Position().Direction();
701         dir2 = Co.Position().XDirection();
702         param1 = apex.Distance(center)*Abs(Tan(angl));
703       }
704       else if (cost < sina ) {
705         typeres = IntAna_Hyperbola;
706         nbint = 2;
707         axex.Normalize();
708
709         deltacenter = sint*sina*sina*distance/(sina*sina - cost*cost);
710         pt1.SetXYZ(center.XYZ() - deltacenter*axex);
711         pt2 = pt1;
712         dir1 = normp;
713         dir2.SetXYZ(axex);
714         param1    = param2 = cost*sina*cosa*distance /(sina*sina-cost*cost);
715         param1bis = param2bis = cost*sina*distance / Sqrt(sina*sina-cost*cost);
716
717       }
718       else {                    // on a alors cost > sina
719         typeres = IntAna_Ellipse;
720         nbint = 1;
721         Standard_Real radius = cost*sina*cosa*distance/(cost*cost-sina*sina);
722         deltacenter = sint*sina*sina*distance/(cost*cost-sina*sina);
723         axex.Normalize();
724         pt1.SetXYZ(center.XYZ() + deltacenter*axex);
725         dir1 = normp;
726         dir2.SetXYZ(axex);
727         param1    = radius;
728         param1bis = cost*sina*distance/ Sqrt(cost*cost - sina*sina);
729       }
730     }
731   }
732   
733   //-- On a du mal a gerer plus loin (Value ProjLib, Params ... )
734   //-- des hyperboles trop bizarres
735   //-- On retourne False -> Traitement par biparametree
736   static Standard_Real EllipseLimit   = 1.0E+9; //OCC513(apo) 1000000
737   static Standard_Real HyperbolaLimit = 2.0E+6; //OCC537(apo) 50000
738   if(typeres==IntAna_Ellipse && nbint>=1) { 
739     if(Abs(param1) > EllipseLimit || Abs(param1bis) > EllipseLimit)  { 
740       done=Standard_False; 
741       return;
742     }
743   }
744   if(typeres==IntAna_Hyperbola && nbint>=2) { 
745     if(Abs(param2) > HyperbolaLimit || Abs(param2bis) > HyperbolaLimit)  { 
746       done = Standard_False; 
747       return;
748     }
749   }
750   if(typeres==IntAna_Hyperbola && nbint>=1) { 
751     if(Abs(param1) > HyperbolaLimit || Abs(param1bis) > HyperbolaLimit)  {
752       done=Standard_False;
753       return;
754     }
755   }
756
757   done = Standard_True;
758 }
759
760 //=======================================================================
761 //function : IntAna_QuadQuadGeo
762 //purpose  : Pln Sphere
763 //=======================================================================
764   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P,
765                                          const gp_Sphere& S)
766 : done(Standard_False),
767   nbint(0),
768   typeres(IntAna_Empty),
769   pt1(0,0,0),
770   pt2(0,0,0),
771   param1(0),
772   param2(0),
773   param1bis(0),
774   param2bis(0),
775   myCommonGen(Standard_False),
776   myPChar(0,0,0)
777 {
778   InitTolerances();
779   Perform(P,S);
780 }
781 //=======================================================================
782 //function : Perform
783 //purpose  : 
784 //=======================================================================
785   void IntAna_QuadQuadGeo::Perform( const gp_Pln& P
786                                    ,const gp_Sphere& S) 
787 {
788   
789   done = Standard_False;
790   Standard_Real A,B,C,D,dist, radius;
791   Standard_Real X,Y,Z;
792
793   nbint = 0;
794 // debug JAG : on met typeres = IntAna_Empty par defaut...
795   typeres = IntAna_Empty;
796   
797   P.Coefficients(A,B,C,D);
798   S.Location().Coord(X,Y,Z);
799   radius = S.Radius();
800   
801   dist = A * X + B * Y + C * Z + D;
802   
803   if (Abs( Abs(dist) - radius) < Epsilon(radius)) {
804     // on a une seule solution : le point projection du centre de la sphere
805     // sur le plan
806     nbint = 1;
807     typeres = IntAna_Point;
808     pt1.SetCoord(X - dist*A, Y - dist*B, Z - dist*C);
809   }
810   else if (Abs(dist) < radius) {
811     // on a un cercle solution
812     nbint = 1;
813     typeres = IntAna_Circle;
814     pt1.SetCoord(X - dist*A, Y - dist*B, Z - dist*C);
815     dir1 = P.Axis().Direction();
816     if(P.Direct()==Standard_False) dir1.Reverse();
817     dir2 = P.Position().XDirection();
818     param1 = Sqrt(radius*radius - dist*dist);
819   }
820   param2bis=0.0; //-- pour eviter param2bis not used .... 
821   done = Standard_True;
822 }
823
824 //=======================================================================
825 //function : IntAna_QuadQuadGeo
826 //purpose  : Cylinder - Cylinder
827 //=======================================================================
828   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl1,
829                                          const gp_Cylinder& Cyl2,
830                                          const Standard_Real Tol) 
831 : done(Standard_False),
832   nbint(0),
833   typeres(IntAna_Empty),
834   pt1(0,0,0),
835   pt2(0,0,0),
836   param1(0),
837   param2(0),
838   param1bis(0),
839   param2bis(0),
840   myCommonGen(Standard_False),
841   myPChar(0,0,0)
842 {
843   InitTolerances();
844   Perform(Cyl1,Cyl2,Tol);
845 }
846 //=======================================================================
847 //function : Perform
848 //purpose  : 
849 //=======================================================================
850   void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl1,
851                      const gp_Cylinder& Cyl2,
852                      const Standard_Real Tol) 
853 {
854   done=Standard_True;
855   //---------------------------- Parallel axes -------------------------
856   AxeOperator A1A2(Cyl1.Axis(),Cyl2.Axis());
857   Standard_Real R1=Cyl1.Radius();
858   Standard_Real R2=Cyl2.Radius();
859   Standard_Real RmR, RmR_Relative;
860   RmR=(R1>R2)? (R1-R2) : (R2-R1);
861   {
862     Standard_Real Rmax, Rmin;
863     Rmax=(R1>R2)? R1 : R2;
864     Rmin=(R1>R2)? R2 : R1;
865     RmR_Relative=RmR/Rmax;
866   }
867
868   Standard_Real DistA1A2=A1A2.Distance();
869   
870   if(A1A2.Parallel()) {
871     if(DistA1A2<=Tol) {
872       if(RmR<=Tol) {
873         typeres=IntAna_Same;
874       }
875       else {
876         typeres=IntAna_Empty;
877       }
878     }
879     else {  //-- DistA1A2 > Tol
880       gp_Pnt P1=Cyl1.Location();
881       gp_Pnt P2t=Cyl2.Location();
882       gp_Pnt P2;
883       //-- P2t is projected on the plane (P1,DirCylX,DirCylY)
884       gp_Dir DirCyl = Cyl1.Position().Direction();
885       Standard_Real ProjP2OnDirCyl1=gp_Vec(DirCyl).Dot(gp_Vec(P1,P2t));
886       
887       P2.SetCoord( P2t.X() - ProjP2OnDirCyl1*DirCyl.X()
888                   ,P2t.Y() - ProjP2OnDirCyl1*DirCyl.Y()
889                   ,P2t.Z() - ProjP2OnDirCyl1*DirCyl.Z());
890       //-- 
891       Standard_Real R1pR2=R1+R2;
892       if(DistA1A2>(R1pR2+Tol)) { 
893         typeres=IntAna_Empty;
894         nbint=0;
895       }
896       else if(DistA1A2>(R1pR2)) {
897         //-- 1 Tangent line -------------------------------------OK
898         typeres=IntAna_Line;
899
900         nbint=1;
901         dir1=DirCyl;
902         Standard_Real R1_R1pR2=R1/R1pR2;
903         pt1.SetCoord( P1.X() + R1_R1pR2 * (P2.X()-P1.X())
904                      ,P1.Y() + R1_R1pR2 * (P2.Y()-P1.Y())
905                      ,P1.Z() + R1_R1pR2 * (P2.Z()-P1.Z()));
906         
907       }
908       else if(DistA1A2>RmR) {
909         //-- 2 lines ---------------------------------------------OK
910         typeres=IntAna_Line;
911         nbint=2;
912         dir1=DirCyl;
913         gp_Vec P1P2(P1,P2);
914         gp_Dir DirA1A2=gp_Dir(P1P2);
915         gp_Dir Ortho_dir1_P1P2 = dir1.Crossed(DirA1A2);
916         dir2=dir1;
917         Standard_Real Alpha=0.5*(R1*R1-R2*R2+DistA1A2*DistA1A2)/(DistA1A2);       
918
919 //      Standard_Real Beta = Sqrt(R1*R1-Alpha*Alpha);
920         Standard_Real anSqrtArg = R1*R1-Alpha*Alpha;
921         Standard_Real Beta = (anSqrtArg > 0.) ? Sqrt(anSqrtArg) : 0.;
922         
923         if((Beta+Beta)<Tol) { 
924           nbint=1;
925           pt1.SetCoord( P1.X() + Alpha*DirA1A2.X()
926                        ,P1.Y() + Alpha*DirA1A2.Y()
927                        ,P1.Z() + Alpha*DirA1A2.Z());
928         }
929         else { 
930           pt1.SetCoord( P1.X() + Alpha*DirA1A2.X() + Beta*Ortho_dir1_P1P2.X()
931                        ,P1.Y() + Alpha*DirA1A2.Y() + Beta*Ortho_dir1_P1P2.Y()
932                        ,P1.Z() + Alpha*DirA1A2.Z() + Beta*Ortho_dir1_P1P2.Z() );
933           
934           pt2.SetCoord( P1.X() + Alpha*DirA1A2.X() - Beta*Ortho_dir1_P1P2.X()
935                        ,P1.Y() + Alpha*DirA1A2.Y() - Beta*Ortho_dir1_P1P2.Y()
936                        ,P1.Z() + Alpha*DirA1A2.Z() - Beta*Ortho_dir1_P1P2.Z());
937         }
938       }
939       else if(DistA1A2>(RmR-Tol)) {
940         //-- 1 Tangent ------------------------------------------OK
941         typeres=IntAna_Line;
942         nbint=1;
943         dir1=DirCyl;
944         Standard_Real R1_RmR=R1/RmR;
945
946         if(R1 < R2) R1_RmR = -R1_RmR;
947
948         pt1.SetCoord( P1.X() + R1_RmR * (P2.X()-P1.X())
949                      ,P1.Y() + R1_RmR * (P2.Y()-P1.Y())
950                      ,P1.Z() + R1_RmR * (P2.Z()-P1.Z()));
951       }
952       else {
953         nbint=0;
954         typeres=IntAna_Empty;
955       }
956     }
957   }
958   else { //-- No Parallel Axis ---------------------------------OK
959     if((RmR_Relative<=myEPSILON_CYLINDER_DELTA_RADIUS) 
960        && (DistA1A2 <= myEPSILON_CYLINDER_DELTA_DISTANCE)) {
961       //-- PI/2 between the two axis   and   Intersection  
962       //-- and identical radius
963       typeres=IntAna_Ellipse;
964       nbint=2;
965       gp_Dir DirCyl1=Cyl1.Position().Direction();
966       gp_Dir DirCyl2=Cyl2.Position().Direction();
967       pt1=pt2=A1A2.PtIntersect();
968       
969       Standard_Real A=DirCyl1.Angle(DirCyl2);
970       Standard_Real B;
971       B=Abs(Sin(0.5*(M_PI-A)));
972       A=Abs(Sin(0.5*A));
973       
974       if(A==0.0 || B==0.0) {
975         typeres=IntAna_Same;
976         return;
977       }
978       
979       
980       gp_Vec dircyl1(DirCyl1);gp_Vec dircyl2(DirCyl2);
981       dir1 = gp_Dir(dircyl1.Added(dircyl2));
982       dir2 = gp_Dir(dircyl1.Subtracted(dircyl2));
983         
984       param2   = Cyl1.Radius() / A;
985       param1   = Cyl1.Radius() / B;
986       param2bis= param1bis = Cyl1.Radius();
987       if(param1 < param1bis) {
988         A=param1; param1=param1bis; param1bis=A;
989       }
990       if(param2 < param2bis) {
991         A=param2; param2=param2bis; param2bis=A;
992       }
993     }
994     else {
995       if(Abs(DistA1A2-Cyl1.Radius()-Cyl2.Radius())<Tol) { 
996         typeres = IntAna_Point;
997         Standard_Real d,p1,p2;
998
999         gp_Dir D1 = Cyl1.Axis().Direction();
1000         gp_Dir D2 = Cyl2.Axis().Direction();
1001         A1A2.Distance(d,p1,p2);
1002         gp_Pnt P = Cyl1.Axis().Location();
1003         gp_Pnt P1(P.X() - p1*D1.X(),
1004                   P.Y() - p1*D1.Y(),
1005                   P.Z() - p1*D1.Z());
1006         P = Cyl2.Axis().Location();
1007         gp_Pnt P2(P.X() - p2*D2.X(),
1008                   P.Y() - p2*D2.Y(),
1009                   P.Z() - p2*D2.Z());
1010         gp_Vec P1P2(P1,P2);
1011         D1=gp_Dir(P1P2);
1012         p1=Cyl1.Radius();
1013         pt1.SetCoord(P1.X() + p1*D1.X(),
1014                      P1.Y() + p1*D1.Y(),
1015                      P1.Z() + p1*D1.Z());
1016         nbint = 1;
1017       }
1018       else {
1019         typeres=IntAna_NoGeometricSolution;
1020       }
1021     }
1022   }
1023 }
1024 //=======================================================================
1025 //function : IntAna_QuadQuadGeo
1026 //purpose  : Cylinder - Cone
1027 //=======================================================================
1028   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
1029                                          const gp_Cone& Con,
1030                                          const Standard_Real Tol) 
1031 : done(Standard_False),
1032   nbint(0),
1033   typeres(IntAna_Empty),
1034   pt1(0,0,0),
1035   pt2(0,0,0),
1036   param1(0),
1037   param2(0),
1038   param1bis(0),
1039   param2bis(0),
1040   myCommonGen(Standard_False),
1041   myPChar(0,0,0)
1042 {
1043   InitTolerances();
1044   Perform(Cyl,Con,Tol);
1045 }
1046 //=======================================================================
1047 //function : Perform
1048 //purpose  : 
1049 //=======================================================================
1050   void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl,
1051                                    const gp_Cone& Con,
1052                                    const Standard_Real ) 
1053 {
1054   done=Standard_True;
1055   AxeOperator A1A2(Cyl.Axis(),Con.Axis());
1056   if(A1A2.Same()) {
1057     gp_Pnt Pt=Con.Apex();
1058     Standard_Real dist=Cyl.Radius()/(Tan(Con.SemiAngle()));
1059     gp_Dir dir=Cyl.Position().Direction();
1060     pt1.SetCoord( Pt.X() + dist*dir.X()
1061                  ,Pt.Y() + dist*dir.Y()
1062                  ,Pt.Z() + dist*dir.Z());
1063     pt2.SetCoord( Pt.X() - dist*dir.X()
1064                  ,Pt.Y() - dist*dir.Y()
1065                  ,Pt.Z() - dist*dir.Z());
1066     dir1=dir2=dir;
1067     param1=param2=Cyl.Radius();
1068     nbint=2;
1069     typeres=IntAna_Circle;
1070
1071   }
1072   else {
1073     typeres=IntAna_NoGeometricSolution;
1074   }
1075 }
1076 //=======================================================================
1077 //function : 
1078 //purpose  : Cylinder - Sphere
1079 //=======================================================================
1080   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
1081                                          const gp_Sphere& Sph,
1082                                          const Standard_Real Tol) 
1083 : done(Standard_False),
1084   nbint(0),
1085   typeres(IntAna_Empty),
1086   pt1(0,0,0),
1087   pt2(0,0,0),
1088   param1(0),
1089   param2(0),
1090   param1bis(0),
1091   param2bis(0),
1092   myCommonGen(Standard_False),
1093   myPChar(0,0,0)
1094 {
1095   InitTolerances();
1096   Perform(Cyl,Sph,Tol);
1097 }
1098 //=======================================================================
1099 //function : Perform
1100 //purpose  : 
1101 //=======================================================================
1102   void IntAna_QuadQuadGeo::Perform( const gp_Cylinder& Cyl
1103                                    ,const gp_Sphere& Sph
1104                                    ,const Standard_Real)
1105 {
1106   done=Standard_True;
1107   gp_Pnt Pt=Sph.Location();
1108   AxeOperator A1A2(Cyl.Axis(),Sph.Position().Axis());
1109   if((A1A2.Intersect()  && Pt.Distance(A1A2.PtIntersect())==0.0 )
1110      || (A1A2.Same()))      {
1111     if(Sph.Radius() < Cyl.Radius()) { 
1112       typeres = IntAna_Empty;
1113     }
1114     else { 
1115       Standard_Real dist=Sqrt( Sph.Radius() * Sph.Radius() - Cyl.Radius() * Cyl.Radius() );
1116       gp_Dir dir=Cyl.Position().Direction();
1117       dir1 = dir2 = dir;
1118       typeres=IntAna_Circle;
1119       pt1.SetCoord( Pt.X() + dist*dir.X()
1120                    ,Pt.Y() + dist*dir.Y()
1121                    ,Pt.Z() + dist*dir.Z());
1122       nbint=1;
1123       param1 = Cyl.Radius();
1124       if(dist>RealEpsilon()) {
1125         pt2.SetCoord( Pt.X() - dist*dir.X()
1126                      ,Pt.Y() - dist*dir.Y()
1127                      ,Pt.Z() - dist*dir.Z());
1128         param2=Cyl.Radius();
1129         nbint=2;
1130       }
1131     }
1132   }
1133   else {
1134     typeres=IntAna_NoGeometricSolution; 
1135   }
1136 }
1137
1138 //=======================================================================
1139 //function : IntAna_QuadQuadGeo
1140 //purpose  : Cone - Cone
1141 //=======================================================================
1142   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cone& Con1,
1143                                          const gp_Cone& Con2,
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   param1(0),
1151   param2(0),
1152   param1bis(0),
1153   param2bis(0),
1154   myCommonGen(Standard_False),
1155   myPChar(0,0,0)
1156 {
1157   InitTolerances();
1158   Perform(Con1,Con2,Tol);
1159 }
1160 //
1161 //=======================================================================
1162 //function : Perform
1163 //purpose  : 
1164 //=======================================================================
1165   void IntAna_QuadQuadGeo::Perform(const gp_Cone& Con1,
1166                                    const gp_Cone& Con2,
1167                                    const Standard_Real Tol) 
1168 {
1169   done=Standard_True;
1170   //
1171   Standard_Real tg1, tg2, aDA1A2, aTol2;
1172   gp_Pnt aPApex1, aPApex2;
1173
1174   Standard_Real TOL_APEX_CONF = 1.e-10;
1175   
1176   //
1177   tg1=Tan(Con1.SemiAngle());
1178   tg2=Tan(Con2.SemiAngle());
1179
1180   if((tg1 * tg2) < 0.) {
1181     tg2 = -tg2;
1182   }
1183   //
1184   aTol2=Tol*Tol;
1185   aPApex1=Con1.Apex();
1186   aPApex2=Con2.Apex();
1187   aDA1A2=aPApex1.SquareDistance(aPApex2);
1188   //
1189   AxeOperator A1A2(Con1.Axis(),Con2.Axis());
1190   //
1191   // 1
1192   if(A1A2.Same()) {
1193     //-- two circles 
1194     Standard_Real x;
1195     gp_Pnt P=Con1.Apex();
1196     gp_Dir D=Con1.Position().Direction();
1197     Standard_Real d=gp_Vec(D).Dot(gp_Vec(P,Con2.Apex()));
1198     
1199     if(Abs(tg1-tg2)>myEPSILON_ANGLE_CONE) { 
1200       if (fabs(d) < TOL_APEX_CONF) {
1201         typeres = IntAna_Point;
1202         nbint = 1;
1203         pt1 = P;
1204         return;
1205       }
1206       x=(d*tg2)/(tg1+tg2);
1207       pt1.SetCoord( P.X() + x*D.X()
1208                    ,P.Y() + x*D.Y()
1209                    ,P.Z() + x*D.Z());
1210       param1=Abs(x*tg1);
1211
1212       x=(d*tg2)/(tg2-tg1);
1213       pt2.SetCoord( P.X() + x*D.X()
1214                    ,P.Y() + x*D.Y()
1215                    ,P.Z() + x*D.Z());
1216       param2=Abs(x*tg1);
1217       dir1 = dir2 = D;
1218       nbint=2;
1219       typeres=IntAna_Circle;
1220     }
1221     else {
1222       if (fabs(d) < TOL_APEX_CONF) { 
1223         typeres=IntAna_Same;
1224       }
1225       else {
1226         typeres=IntAna_Circle;
1227         nbint=1;
1228         x=d*0.5;
1229         pt1.SetCoord( P.X() + x*D.X()
1230                      ,P.Y() + x*D.Y()
1231                      ,P.Z() + x*D.Z());
1232         param1 = Abs(x * tg1);
1233         dir1 = D;
1234       }
1235     }
1236   } //-- fin A1A2.Same
1237   // 2
1238   else if((Abs(tg1-tg2)<myEPSILON_ANGLE_CONE) && (A1A2.Parallel())) {
1239     //-- voir AnVer12mai98
1240     Standard_Real DistA1A2=A1A2.Distance();
1241     gp_Dir DA1=Con1.Position().Direction();
1242     gp_Vec O1O2(Con1.Apex(),Con2.Apex());
1243     Standard_Real O1O2_DA1=gp_Vec(DA1).Dot(O1O2);
1244     
1245     gp_Vec O1_Proj_A2(O1O2.X()-O1O2_DA1*DA1.X(),
1246                       O1O2.Y()-O1O2_DA1*DA1.Y(),
1247                       O1O2.Z()-O1O2_DA1*DA1.Z());
1248     gp_Dir DB1=gp_Dir(O1_Proj_A2);
1249     
1250     Standard_Real yO1O2=O1O2.Dot(gp_Vec(DA1));
1251     Standard_Real ABSTG1 = Abs(tg1);
1252     Standard_Real X2 = (DistA1A2/ABSTG1 - yO1O2)*0.5;
1253     Standard_Real X1 = X2+yO1O2;
1254     
1255     gp_Pnt P1(Con1.Apex().X() + X1*( DA1.X() + ABSTG1*DB1.X()),
1256               Con1.Apex().Y() + X1*( DA1.Y() + ABSTG1*DB1.Y()), 
1257               Con1.Apex().Z() + X1*( DA1.Z() + ABSTG1*DB1.Z()));
1258
1259     gp_Pnt MO1O2(0.5*(Con1.Apex().X()+Con2.Apex().X()),
1260                  0.5*(Con1.Apex().Y()+Con2.Apex().Y()),
1261                  0.5*(Con1.Apex().Z()+Con2.Apex().Z()));
1262     gp_Vec P1MO1O2(P1,MO1O2);
1263     
1264     gp_Dir DA1_X_DB1=DA1.Crossed(DB1);
1265     gp_Dir OrthoPln =  DA1_X_DB1.Crossed(gp_Dir(P1MO1O2));
1266     
1267     IntAna_QuadQuadGeo INTER_QUAD_PLN(gp_Pln(P1,OrthoPln),Con1,Tol,Tol);
1268       if(INTER_QUAD_PLN.IsDone()) {
1269       switch(INTER_QUAD_PLN.TypeInter()) {
1270       case IntAna_Ellipse:      {
1271         typeres=IntAna_Ellipse;
1272         gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
1273         pt1 = E.Location();
1274         dir1 = E.Position().Direction();
1275         dir2 = E.Position().XDirection();
1276         param1 = E.MajorRadius();
1277         param1bis = E.MinorRadius();
1278         nbint = 1;
1279         break; 
1280       }
1281       case IntAna_Circle: {
1282         typeres=IntAna_Circle;
1283         gp_Circ C=INTER_QUAD_PLN.Circle(1);
1284         pt1 = C.Location();
1285         dir1 = C.Position().XDirection();
1286         dir2 = C.Position().YDirection();
1287         param1 = C.Radius();
1288         nbint = 1;
1289         break;
1290       }
1291       case IntAna_Hyperbola: {
1292         typeres=IntAna_Hyperbola;
1293         gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
1294         pt1 = pt2 = H.Location();
1295         dir1 = H.Position().Direction();
1296         dir2 = H.Position().XDirection();
1297         param1 = param2 = H.MajorRadius();
1298         param1bis = param2bis = H.MinorRadius();
1299         nbint = 2;
1300         break;
1301       }
1302       case IntAna_Line: {
1303         typeres=IntAna_Line;
1304         gp_Lin H=INTER_QUAD_PLN.Line(1);
1305         pt1 = pt2 = H.Location();
1306         dir1 = dir2 = H.Position().Direction();
1307         param1 = param2 = 0.0;
1308         param1bis = param2bis = 0.0;
1309         nbint = 2;
1310         break;
1311       }
1312       default:
1313         typeres=IntAna_NoGeometricSolution; 
1314       }
1315     }
1316   }// else if((Abs(tg1-tg2)<EPSILON_ANGLE_CONE) && (A1A2.Parallel()))
1317   // 3
1318   else if (aDA1A2<aTol2) {
1319     //
1320     // When apices are coinsided there can be 3 possible cases
1321     // 3.1 - empty solution (iRet=0)
1322     // 3.2 - one line  when cone1 touches cone2 (iRet=1)
1323     // 3.3 - two lines when cone1 intersects cone2 (iRet=2)
1324     //
1325     Standard_Integer iRet;
1326     Standard_Real aGamma, aBeta1, aBeta2;
1327     Standard_Real aD1, aR1, aTgBeta1, aTgBeta2, aHalfPI;
1328     Standard_Real aCosGamma, aSinGamma, aDx, aR2, aRD2, aD2;
1329     gp_Pnt2d aP0, aPA1, aP1, aPA2;
1330     gp_Vec2d aVAx2;
1331     gp_Ax1 aAx1, aAx2;
1332     //
1333     // Preliminary analysis. Determination of iRet
1334     //
1335     iRet=0;
1336     aHalfPI=0.5*M_PI;
1337     aD1=1.;
1338     aPA1.SetCoord(aD1, 0.);
1339     aP0.SetCoord(0., 0.);
1340     //
1341     aAx1=Con1.Axis();
1342     aAx2=Con2.Axis();
1343     aGamma=aAx1.Angle(aAx2);
1344     if (aGamma>aHalfPI){
1345       aGamma=M_PI-aGamma;
1346     }
1347     aCosGamma=Cos(aGamma);
1348     aSinGamma=Sin(aGamma);
1349     //
1350     aBeta1=Con1.SemiAngle();
1351     aTgBeta1=Tan(aBeta1);
1352     aTgBeta1=Abs(aTgBeta1);
1353     //
1354     aBeta2=Con2.SemiAngle();
1355     aTgBeta2=Tan(aBeta2);
1356     aTgBeta2=Abs(aTgBeta2);
1357     //
1358     aR1=aD1*aTgBeta1;
1359     aP1.SetCoord(aD1, aR1);
1360     //
1361     // PA2
1362     aVAx2.SetCoord(aCosGamma, aSinGamma);
1363     gp_Dir2d aDAx2(aVAx2);
1364     gp_Lin2d aLAx2(aP0, aDAx2);
1365     //
1366     gp_Vec2d aV(aP0, aP1);
1367     aDx=aVAx2.Dot(aV);
1368     aPA2=aP0.Translated(aDx*aDAx2);
1369     //
1370     // aR2
1371     aDx=aPA2.Distance(aP0);
1372     aR2=aDx*aTgBeta2;
1373     //
1374     // aRD2
1375     aRD2=aPA2.Distance(aP1);
1376     //
1377     if (aRD2>(aR2+Tol)) {
1378       iRet=0;
1379       typeres=IntAna_Empty; //nothing
1380       return;
1381     }
1382     //
1383     iRet=1; //touch case => 1 line
1384     if (aRD2<(aR2-Tol)) {
1385       iRet=2;//intersection => couple of lines
1386     }
1387     //
1388     // Finding the solution in 3D
1389     //
1390     Standard_Real aDa;
1391     gp_Pnt aQApex1, aQA1, aQA2, aQX, aQX1, aQX2;
1392     gp_Dir aD3Ax1, aD3Ax2;
1393     gp_Lin aLin;
1394     IntAna_QuadQuadGeo aIntr;
1395     //
1396     aQApex1=Con1.Apex();
1397     aD3Ax1=aAx1.Direction(); 
1398     aQA1.SetCoord(aQApex1.X()+aD1*aD3Ax1.X(),
1399                   aQApex1.Y()+aD1*aD3Ax1.Y(),
1400                   aQApex1.Z()+aD1*aD3Ax1.Z());
1401     //
1402     aDx=aD3Ax1.Dot(aAx2.Direction());
1403     if (aDx<0.) {
1404       aAx2.Reverse();
1405     }
1406     aD3Ax2=aAx2.Direction();
1407     //
1408     aD2=aD1*sqrt((1.+aTgBeta1*aTgBeta1)/(1.+aTgBeta2*aTgBeta2));
1409     //
1410     aQA2.SetCoord(aQApex1.X()+aD2*aD3Ax2.X(),
1411                   aQApex1.Y()+aD2*aD3Ax2.Y(),
1412                   aQApex1.Z()+aD2*aD3Ax2.Z());
1413     //
1414     gp_Pln aPln1(aQA1, aD3Ax1);
1415     gp_Pln aPln2(aQA2, aD3Ax2);
1416     //
1417     aIntr.Perform(aPln1, aPln2, Tol, Tol);
1418     if (!aIntr.IsDone()) {
1419       iRet=-1; // just in case. it must not be so
1420       typeres=IntAna_NoGeometricSolution; 
1421       return;
1422     }
1423     //
1424     aLin=aIntr.Line(1);
1425     const gp_Dir& aDLin=aLin.Direction();
1426     gp_Vec aVLin(aDLin);
1427     gp_Pnt aOrig=aLin.Location();
1428     gp_Vec aVr(aQA1, aOrig);
1429     aDx=aVLin.Dot(aVr);
1430     aQX=aOrig.Translated(aDx*aVLin);
1431     //
1432     // Final part
1433     //
1434     typeres=IntAna_Line; 
1435     //
1436     param1=0.;
1437     param2 =0.;
1438     param1bis=0.;
1439     param2bis=0.;
1440     //
1441     if (iRet==1) {
1442       // one line
1443       nbint=1;
1444       pt1=aQApex1;
1445       gp_Vec aVX(aQApex1, aQX);
1446       dir1=gp_Dir(aVX);
1447     }
1448     
1449     else {//iRet=2 
1450       // two lines
1451       nbint=2;
1452       aDa=aQA1.Distance(aQX);
1453       aDx=sqrt(aR1*aR1-aDa*aDa);
1454       aQX1=aQX.Translated(aDx*aVLin);
1455       aQX2=aQX.Translated(-aDx*aVLin);
1456       //
1457       pt1=aQApex1;
1458       pt2=aQApex1;
1459       gp_Vec aVX1(aQApex1, aQX1);
1460       dir1=gp_Dir(aVX1);
1461       gp_Vec aVX2(aQApex1, aQX2);
1462       dir2=gp_Dir(aVX2);
1463     }
1464   } //else if (aDA1A2<aTol2) {
1465   //Case when cones have common generatrix
1466   else if(A1A2.Intersect()) {
1467     //Check if apex of one cone belongs another one
1468     Standard_Real u, v, tol2 = Tol*Tol;
1469     ElSLib::Parameters(Con2, aPApex1, u, v);
1470     gp_Pnt p = ElSLib::Value(u, v, Con2);
1471     if(aPApex1.SquareDistance(p) > tol2) {
1472       typeres=IntAna_NoGeometricSolution; 
1473       return;
1474     }
1475     //
1476     ElSLib::Parameters(Con1, aPApex2, u, v);
1477     p = ElSLib::Value(u, v, Con1);
1478     if(aPApex2.SquareDistance(p) > tol2) {
1479       typeres=IntAna_NoGeometricSolution; 
1480       return;
1481     }
1482
1483     //Cones have a common generatrix passing through apexes
1484     myCommonGen = Standard_True;
1485
1486     //common generatrix of cones
1487     gp_Lin aGen(aPApex1, gp_Dir(gp_Vec(aPApex1, aPApex2)));
1488
1489     //Intersection point of axes
1490     gp_Pnt aPAxeInt = A1A2.PtIntersect();
1491
1492     //Characteristic point of intersection curve
1493     u = ElCLib::Parameter(aGen, aPAxeInt);
1494     myPChar = ElCLib::Value(u, aGen);
1495
1496
1497     //Other generatrixes of cones laying in maximal plane
1498     gp_Lin aGen1 = aGen.Rotated(Con1.Axis(), M_PI); 
1499     gp_Lin aGen2 = aGen.Rotated(Con2.Axis(), M_PI); 
1500     //
1501     //Intersection point of generatrixes
1502     gp_Dir aN; //solution plane normal
1503     gp_Dir aD1 = aGen1.Direction();
1504
1505     gp_Dir aD2(aD1.Crossed(aGen.Direction()));
1506
1507     if(aD1.IsParallel(aGen2.Direction(), Precision::Angular())) {
1508       aN = aD1.Crossed(aD2);
1509     }
1510     else if(aGen1.SquareDistance(aGen2) > tol2) {
1511       //Something wrong ???
1512       typeres=IntAna_NoGeometricSolution; 
1513       return;
1514     }
1515     else {
1516       gp_Dir D1 = aGen1.Position().Direction();
1517       gp_Dir D2 = aGen2.Position().Direction();
1518       gp_Pnt O1 = aGen1.Location();
1519       gp_Pnt O2 = aGen2.Location();
1520       Standard_Real D1DotD2 = D1.Dot(D2);
1521       Standard_Real aSin = 1.-D1DotD2*D1DotD2;
1522       gp_Vec O1O2 (O1,O2);
1523       Standard_Real U2 = (D1.XYZ()*(O1O2.Dot(D1))-(O1O2.XYZ())).Dot(D2.XYZ());
1524       U2 /= aSin;
1525       gp_Pnt aPGint(ElCLib::Value(U2, aGen2));
1526     
1527       aD1 = gp_Dir(gp_Vec(aPGint, myPChar));
1528       aN = aD1.Crossed(aD2);
1529     }
1530     //Plane that must contain intersection curves
1531     gp_Pln anIntPln(myPChar, aN);
1532
1533     IntAna_QuadQuadGeo INTER_QUAD_PLN(anIntPln,Con1,Tol,Tol);
1534
1535       if(INTER_QUAD_PLN.IsDone()) {
1536       switch(INTER_QUAD_PLN.TypeInter()) {
1537       case IntAna_Ellipse:      {
1538         typeres=IntAna_Ellipse;
1539         gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
1540         pt1 = E.Location();
1541         dir1 = E.Position().Direction();
1542         dir2 = E.Position().XDirection();
1543         param1 = E.MajorRadius();
1544         param1bis = E.MinorRadius();
1545         nbint = 1;
1546         break; 
1547       }
1548       case IntAna_Circle: {
1549         typeres=IntAna_Circle;
1550         gp_Circ C=INTER_QUAD_PLN.Circle(1);
1551         pt1 = C.Location();
1552         dir1 = C.Position().XDirection();
1553         dir2 = C.Position().YDirection();
1554         param1 = C.Radius();
1555         nbint = 1;
1556         break;
1557       }
1558       case IntAna_Parabola: {
1559         typeres=IntAna_Parabola;
1560         gp_Parab Prb=INTER_QUAD_PLN.Parabola(1);
1561         pt1 = Prb.Location();
1562         dir1 = Prb.Position().Direction();
1563         dir2 = Prb.Position().XDirection();
1564         param1 = Prb.Focal();
1565         nbint = 1;
1566         break;
1567       }
1568       case IntAna_Hyperbola: {
1569         typeres=IntAna_Hyperbola;
1570         gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
1571         pt1 = pt2 = H.Location();
1572         dir1 = H.Position().Direction();
1573         dir2 = H.Position().XDirection();
1574         param1 = param2 = H.MajorRadius();
1575         param1bis = param2bis = H.MinorRadius();
1576         nbint = 2;
1577         break;
1578       }
1579       default:
1580         typeres=IntAna_NoGeometricSolution; 
1581       }
1582     }
1583   }
1584   
1585   else {
1586     typeres=IntAna_NoGeometricSolution; 
1587   }
1588 }
1589 //=======================================================================
1590 //function : IntAna_QuadQuadGeo
1591 //purpose  : Sphere - Cone
1592 //=======================================================================
1593   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Sphere& Sph,
1594                                          const gp_Cone& Con,
1595                                          const Standard_Real Tol) 
1596 : done(Standard_False),
1597   nbint(0),
1598   typeres(IntAna_Empty),
1599   pt1(0,0,0),
1600   pt2(0,0,0),
1601   param1(0),
1602   param2(0),
1603   param1bis(0),
1604   param2bis(0),
1605   myCommonGen(Standard_False),
1606   myPChar(0,0,0)
1607 {
1608   InitTolerances();
1609   Perform(Sph,Con,Tol);
1610 }
1611 //=======================================================================
1612 //function : Perform
1613 //purpose  : 
1614 //=======================================================================
1615   void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph,
1616                                    const gp_Cone& Con,
1617                                    const Standard_Real)
1618 {
1619   
1620   //
1621   done=Standard_True;
1622   //
1623   AxeOperator A1A2(Con.Axis(),Sph.Position().Axis());
1624   gp_Pnt Pt=Sph.Location();
1625   //
1626   if((A1A2.Intersect() && (Pt.Distance(A1A2.PtIntersect())==0.0))
1627      || A1A2.Same()) {
1628     gp_Pnt ConApex= Con.Apex();
1629     Standard_Real dApexSphCenter=Pt.Distance(ConApex); 
1630     gp_Dir ConDir;
1631     if(dApexSphCenter>RealEpsilon()) { 
1632       ConDir = gp_Dir(gp_Vec(ConApex,Pt));
1633     }
1634     else { 
1635       ConDir = Con.Position().Direction();
1636     }
1637     
1638     Standard_Real Rad=Sph.Radius();
1639     Standard_Real tga=Tan(Con.SemiAngle());
1640
1641
1642     //-- 2 circles
1643     //-- x: Roots of    (x**2 + y**2 = Rad**2)
1644     //--                tga = y / (x+dApexSphCenter)
1645     Standard_Real tgatga = tga * tga;
1646     math_DirectPolynomialRoots Eq( 1.0+tgatga
1647                                   ,2.0*tgatga*dApexSphCenter
1648                                   ,-Rad*Rad + dApexSphCenter*dApexSphCenter*tgatga);
1649     if(Eq.IsDone()) {
1650       Standard_Integer nbsol=Eq.NbSolutions();
1651       if(nbsol==0) {
1652         typeres=IntAna_Empty;
1653       }
1654       else { 
1655         typeres=IntAna_Circle;
1656         if(nbsol>=1) {
1657           Standard_Real x = Eq.Value(1);
1658           Standard_Real dApexSphCenterpx = dApexSphCenter+x;
1659           nbint=1;
1660           pt1.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
1661                        ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
1662                        ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
1663           param1 = tga * dApexSphCenterpx;
1664           param1 = Abs(param1);
1665           dir1 = ConDir;
1666           if(param1<=myEPSILON_MINI_CIRCLE_RADIUS) {
1667             typeres=IntAna_PointAndCircle;
1668             param1=0.0;
1669           }
1670         }
1671         if(nbsol>=2) {
1672           Standard_Real x=Eq.Value(2);
1673           Standard_Real dApexSphCenterpx = dApexSphCenter+x;
1674           nbint=2;
1675           pt2.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
1676                        ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
1677                        ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
1678           param2 = tga * dApexSphCenterpx;
1679           param2 = Abs(param2);
1680           dir2=ConDir;
1681           if(param2<=myEPSILON_MINI_CIRCLE_RADIUS) {
1682             typeres=IntAna_PointAndCircle;
1683             param2=0.0;
1684           }
1685         }
1686       }
1687     }
1688     else {
1689       done=Standard_False;
1690     }
1691   }
1692   else {
1693     typeres=IntAna_NoGeometricSolution; 
1694   }
1695 }
1696
1697 //=======================================================================
1698 //function : IntAna_QuadQuadGeo
1699 //purpose  : Sphere - Sphere
1700 //=======================================================================
1701   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo( const gp_Sphere& Sph1
1702                                          ,const gp_Sphere& Sph2
1703                                          ,const Standard_Real Tol) 
1704 : done(Standard_False),
1705   nbint(0),
1706   typeres(IntAna_Empty),
1707   pt1(0,0,0),
1708   pt2(0,0,0),
1709   param1(0),
1710   param2(0),
1711   param1bis(0),
1712   param2bis(0),
1713   myCommonGen(Standard_False),
1714   myPChar(0,0,0)
1715 {
1716   InitTolerances();
1717   Perform(Sph1,Sph2,Tol);
1718 }
1719 //=======================================================================
1720 //function : Perform
1721 //purpose  : 
1722 //=======================================================================
1723   void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph1,
1724                                    const gp_Sphere& Sph2,
1725                                    const Standard_Real Tol)   
1726 {
1727   done=Standard_True;
1728   gp_Pnt O1=Sph1.Location();
1729   gp_Pnt O2=Sph2.Location();
1730   Standard_Real dO1O2=O1.Distance(O2);
1731   Standard_Real R1=Sph1.Radius();
1732   Standard_Real R2=Sph2.Radius();
1733   Standard_Real Rmin,Rmax;
1734   typeres=IntAna_Empty;
1735   param2bis=0.0; //-- pour eviter param2bis not used .... 
1736
1737   if(R1>R2) { Rmin=R2; Rmax=R1; } else { Rmin=R1; Rmax=R2; }
1738   
1739   if(dO1O2<=Tol && (Abs(R1-R2) <= Tol)) {
1740     typeres = IntAna_Same;
1741   }
1742   else { 
1743     if(dO1O2<=Tol) { return; } 
1744     gp_Dir Dir=gp_Dir(gp_Vec(O1,O2));
1745     Standard_Real t = Rmax - dO1O2 - Rmin;
1746
1747     //----------------------------------------------------------------------
1748     //--        |----------------- R1 --------------------|
1749     //--        |----dO1O2-----|-----------R2----------|
1750     //--                                            --->--<-- t
1751     //--
1752     //--        |------ R1 ------|---------dO1O2----------|
1753     //--     |-------------------R2-----------------------|
1754     //--  --->--<-- t
1755     //----------------------------------------------------------------------
1756     if(t >= 0.0  && t <=Tol) { 
1757       typeres = IntAna_Point;
1758       nbint = 1;
1759       Standard_Real t2;
1760       if(R1==Rmax) t2=(R1 + (R2 + dO1O2)) * 0.5;
1761       else         t2=(-R1+(dO1O2-R2))*0.5;
1762         
1763       pt1.SetCoord( O1.X() + t2*Dir.X()
1764                    ,O1.Y() + t2*Dir.Y()
1765                    ,O1.Z() + t2*Dir.Z());
1766     }
1767     else  {
1768       //-----------------------------------------------------------------
1769       //--        |----------------- dO1O2 --------------------|
1770       //--        |----R1-----|-----------R2----------|-Tol-|
1771       //--                                            
1772       //--        |----------------- Rmax --------------------|
1773       //--        |----Rmin----|-------dO1O2-------|-Tol-|
1774       //--                                            
1775       //-----------------------------------------------------------------
1776       if((dO1O2 > (R1+R2+Tol)) || (Rmax > (dO1O2+Rmin+Tol))) {
1777         typeres=IntAna_Empty;
1778       }
1779       else {
1780         //---------------------------------------------------------------
1781         //--     
1782         //--
1783         //---------------------------------------------------------------
1784         Standard_Real Alpha=0.5*(R1*R1-R2*R2+dO1O2*dO1O2)/(dO1O2);       
1785         Standard_Real Beta = R1*R1-Alpha*Alpha;
1786         Beta = (Beta>0.0)? Sqrt(Beta) : 0.0;
1787         
1788         if(Beta<= myEPSILON_MINI_CIRCLE_RADIUS) { 
1789           typeres = IntAna_Point;
1790           Alpha = (R1 + (dO1O2 - R2)) * 0.5;
1791         }
1792         else { 
1793           typeres = IntAna_Circle;
1794           dir1 = Dir;
1795           param1 = Beta;
1796         }         
1797         pt1.SetCoord( O1.X() + Alpha*Dir.X()
1798                      ,O1.Y() + Alpha*Dir.Y()
1799                      ,O1.Z() + Alpha*Dir.Z());
1800         
1801         nbint=1;
1802       }
1803     }
1804   }
1805 }
1806 //=======================================================================
1807 //function : Point
1808 //purpose  : Returns a Point
1809 //=======================================================================
1810   gp_Pnt IntAna_QuadQuadGeo::Point(const Standard_Integer n) const 
1811 {
1812   if(!done)          {    StdFail_NotDone::Raise();        }
1813   if(n>nbint || n<1) {    Standard_DomainError::Raise();   }
1814   if(typeres==IntAna_PointAndCircle) {
1815     if(n!=1) { Standard_DomainError::Raise();  }
1816     if(param1==0.0) return(pt1);
1817     return(pt2);
1818   }
1819   else if(typeres==IntAna_Point) {
1820     if(n==1) return(pt1);
1821     return(pt2);
1822   }
1823
1824   // WNT (what can you expect from MicroSoft ?)
1825   return gp_Pnt(0,0,0);
1826 }
1827 //=======================================================================
1828 //function : Line
1829 //purpose  : Returns a Line
1830 //=======================================================================
1831   gp_Lin IntAna_QuadQuadGeo::Line(const Standard_Integer n) const 
1832 {
1833   if(!done)        {   StdFail_NotDone::Raise();   }
1834   if((n>nbint) || (n<1) || (typeres!=IntAna_Line)) {
1835     Standard_DomainError::Raise();
1836     }
1837   if(n==1) {  return(gp_Lin(pt1,dir1));   }
1838   else {      return(gp_Lin(pt2,dir2));   }
1839 }
1840 //=======================================================================
1841 //function : Circle
1842 //purpose  : Returns a Circle
1843 //=======================================================================
1844   gp_Circ IntAna_QuadQuadGeo::Circle(const Standard_Integer n) const 
1845 {
1846   if(!done) {    StdFail_NotDone::Raise();     }
1847   if(typeres==IntAna_PointAndCircle) {
1848     if(n!=1) { Standard_DomainError::Raise();  }
1849     if(param2==0.0) return(gp_Circ(DirToAx2(pt1,dir1),param1));
1850     return(gp_Circ(DirToAx2(pt2,dir2),param2));
1851   }
1852   else if((n>nbint) || (n<1) || (typeres!=IntAna_Circle)) {
1853     Standard_DomainError::Raise();
1854     }
1855   if(n==1) { return(gp_Circ(DirToAx2(pt1,dir1),param1));   }
1856   else {     return(gp_Circ(DirToAx2(pt2,dir2),param2));   }
1857 }
1858
1859 //=======================================================================
1860 //function : Ellipse
1861 //purpose  : Returns a Elips  
1862 //=======================================================================
1863   gp_Elips IntAna_QuadQuadGeo::Ellipse(const Standard_Integer n) const
1864 {
1865   if(!done) {     StdFail_NotDone::Raise();     }
1866   if((n>nbint) || (n<1) || (typeres!=IntAna_Ellipse)) {
1867     Standard_DomainError::Raise();
1868   }
1869
1870   if(n==1) {
1871     Standard_Real R1=param1, R2=param1bis, aTmp;
1872     if (R1<R2) {
1873       aTmp=R1; R1=R2; R2=aTmp;
1874     }
1875     gp_Ax2 anAx2(pt1, dir1 ,dir2);
1876     gp_Elips anElips (anAx2, R1, R2);
1877     return anElips;
1878   }
1879   else {
1880     Standard_Real R1=param2, R2=param2bis, aTmp;
1881     if (R1<R2) {
1882       aTmp=R1; R1=R2; R2=aTmp;
1883     }
1884     gp_Ax2 anAx2(pt2, dir2 ,dir1);
1885     gp_Elips anElips (anAx2, R1, R2);
1886     return anElips;
1887   }
1888 }
1889 //=======================================================================
1890 //function : Parabola
1891 //purpose  : Returns a Parabola 
1892 //=======================================================================
1893   gp_Parab IntAna_QuadQuadGeo::Parabola(const Standard_Integer n) const 
1894 {
1895   if(!done) {
1896     StdFail_NotDone::Raise();
1897     }
1898   if (typeres!=IntAna_Parabola) {
1899     Standard_DomainError::Raise();
1900   }
1901   if((n>nbint) || (n!=1)) {
1902     Standard_OutOfRange::Raise();
1903   }
1904   return(gp_Parab(gp_Ax2( pt1
1905                          ,dir1
1906                          ,dir2)
1907                   ,param1));
1908 }
1909 //=======================================================================
1910 //function : Hyperbola
1911 //purpose  : Returns a Hyperbola  
1912 //=======================================================================
1913   gp_Hypr IntAna_QuadQuadGeo::Hyperbola(const Standard_Integer n) const 
1914 {
1915   if(!done) {
1916     StdFail_NotDone::Raise();
1917     }
1918   if((n>nbint) || (n<1) || (typeres!=IntAna_Hyperbola)) {
1919     Standard_DomainError::Raise();
1920     }
1921   if(n==1) {
1922     return(gp_Hypr(gp_Ax2( pt1
1923                           ,dir1
1924                           ,dir2)
1925                    ,param1,param1bis));
1926   }
1927   else {
1928     return(gp_Hypr(gp_Ax2( pt2
1929                           ,dir1
1930                           ,dir2.Reversed())
1931                    ,param2,param2bis));
1932   }
1933 }
1934 //=======================================================================
1935 //function : HasCommonGen
1936 //purpose  : 
1937 //=======================================================================
1938 Standard_Boolean IntAna_QuadQuadGeo::HasCommonGen() const
1939 {
1940   return myCommonGen;
1941 }
1942 //=======================================================================
1943 //function : PChar
1944 //purpose  : 
1945 //=======================================================================
1946 const gp_Pnt& IntAna_QuadQuadGeo::PChar() const
1947 {
1948   return myPChar;
1949 }
1950 //=======================================================================
1951 //function : RefineDir
1952 //purpose  : 
1953 //=======================================================================
1954 void RefineDir(gp_Dir& aDir)
1955 {
1956   Standard_Integer k, m, n;
1957   Standard_Real aC[3];
1958   //
1959   aDir.Coord(aC[0], aC[1], aC[2]);
1960   //
1961   m=0;
1962   n=0;
1963   for (k=0; k<3; ++k) {
1964     if (aC[k]==1. || aC[k]==-1.) {
1965       ++m;
1966     }
1967     else if (aC[k]!=0.) {
1968       ++n;
1969     }
1970   }
1971   //
1972   if (m && n) {
1973     Standard_Real aEps, aR1, aR2, aNum;
1974     //
1975     aEps=RealEpsilon();
1976     aR1=1.-aEps;
1977     aR2=1.+aEps;
1978     //
1979     for (k=0; k<3; ++k) {
1980       m=(aC[k]>0.);
1981       aNum=(m)? aC[k] : -aC[k];
1982       if (aNum>aR1 && aNum<aR2) {
1983         if (m) {
1984           aC[k]=1.;
1985         }         
1986         else {
1987           aC[k]=-1.;
1988         }
1989         //
1990         aC[(k+1)%3]=0.;
1991         aC[(k+2)%3]=0.;
1992         break;
1993       }
1994     }
1995     aDir.SetCoord(aC[0], aC[1], aC[2]);
1996   }
1997 }
1998
1999
2000