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