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