src/IntAna/IntAna_QuadQuadGeo.cxx

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