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