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