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