src/IntAna/IntAna_QuadQuadGeo.cxx

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